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ABSTRACT 



The nontrivial topological structure of the QCD gauge vacuum generates a CP 
breaking term in the QCD Lagrangian. However, measurements of the neutron elec- 
tric dipole moment have demonstrated that the term's coefficient is unnaturally small, 
a dilemma known as the strong CP problem. A massless up quark has long been seen 
as a potential solution, as the term could then be absorbed through the resulting free- 
dom to perform arbitrary chiral rotations on the up quark field. 

Through the light-quark-mass ratio niu/md, leading order Chiral Perturbation 
Theory appears to rule this scenario out. However, the Kaplan-Manohar ambigu- 
ity demonstrates that certain strong next-to-leading order corrections are indistin- 
guishable from the effects of an up quark mass. Only a direct calculation of the 
Gasser-Leutwyler coefficient combination 2L8 — L5 can resolve the issue. 

New theoretical insights into partial quenched Chiral Perturbation Theory have 
revealed that a calculation of the low-energy constants of the partially quenched chiral 
Lagrangian is equivalent to a determination of the physical Gasser-Leutwyler coeffi- 
cients. The coefficient combination in question is directly accessible through the pion 
mass's dependence on the valence quark mass, a dependence ripe for determination 
via Lattice Quantum Chromodynamics. 
ii 



We carry out such a partially quenched lattice calculation using Nf — 3 staggered 
fermions and the recently developed smearing technique known as hypercubic block- 
ing. Through the use of several ensembles, we make a quantitative assessment of our 
systematic error. We find 2Ls — L5 = (0.22 ± 0.14) x 10~^, which corresponds to 
a light-quark-mass ratio of rriu/md — 0.408 ± 0.035. Thus, our study rules out the 
massless-up-quark solution to the strong CP problem. 

This is the first calculation of its type to use a physical number of light quarks, 
Nf — 3, and the first determination of 2Ls — I/5 to include a comprehensive study of 
statistical error. 
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CHAPTER 1 



INTRODUCTION 



Since its inception, the low-energy dynamics of Quantum Chromodynamics (QCD) 
have been poorly understood. To date, the most successful framework for building 
an understanding of low-energy QCD has been Chiral Perturbation Theory, with the 
current large uncertainty in its NLO coefficients quantifying our lackluster under- 
standing. These low-energy constants, the Gasser-Leutwyler coefficients, have error 
bars which range from 10% to 160% of their value. For many of these coefficients, 
their uncertainty has never been reduced below the magnitude determined at the first 
instance of their calculation some twenty years ago. 

As an effective field theory, light-meson Chiral Perturbation Theory collects the 
interactions of QCD into a finite number of meson vertices. The complex low-energy 
dynamics of QCD are boiled down into the coefficients of these vertices, the Gasser- 
Leutwyler coefficients. Thus, while various theoretical and phenomenological methods 
have been used to estimate their values, the most direct determination of the Gasser- 
Leutwyler coefficients would be a calculation of these vertices' strength using the 
fundamental theory, QCD. 

While perturbation theory clearly fails in this regard, as the strong coupling at 
these energy scales is of order one, lattice techniques have proven successful. In fact, 
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this may be a case in which lattice field theory can provide the physics community 
with the best results available. 

We present here a lattice calculation of a single important combination of the 
Gasser-Leutwyler coefficients, 2L^ — L5. The motivation behind this combination 
choice is twofold. First, it is the combination whose calculation on the lattice is 
most straightforward. Additionally, this combination provides insight into the NLO 
corrections to the light-quark-mass ratio m„/mrf. Thus, determining 2L^ — L5, even 
with moderate accuracy, allows one to lay to rest the possibility of a massless up 
quark. 

While we calculate definitively only a single combination of the Gasser-Leutwyler 
coefficients, this study is only the first step in a larger effort by the lattice community 
to generate results for all accessible Gasser-Leutwyler coefficients. 

Elements of this study have seen earlier publication A similar investigation, 
which uses an unphysical number of light quarks and less sophisticated analysis tech- 
niques, can be found in [^. 

In Chapter H, after a brief overview of QCD, the strong CP problem and its 
potential solutions are presented, including the directly relevant solution involving a 
massless up quark. 

In Chapter ^we introduce Chiral Perturbation Theory up to NLO, focusing on the 
insight it imparts into the light-quark-mass ratio m^j/m^^. The relationship between 
the quark-mass ratio and the light mesons is explored, and the importance of the 
Gasser-Leutwyler coefficient combination 2L^ — L5 in that relationship is presented. 
An overview of past phenomenological and theoretical estimates for the coefficients 
L5 and Lg is given, including a discussion of the Kaplan-Manohar ambiguity, which 
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makes an experimental determination of 2^8 — L5 impossible and theoretical estimates 
challenging. 

In Chapter ^ the basics of Lattice Quantum Field Theory are presented. 

In Chapter ^ we build on those basics, introducing Lattice Quantum Chromody- 
namics. The lattice techniques used in this study are covered, including staggered 
fermions, the conjugate gradient method, the R algorithm, hypercubic blocking, and 
the Sommer scale. We present in explicit detail the methods used to extract the pion 
mass and decay constant from staggered bilinear correlators. 

In Chapter |^ we extend the concepts of Chiral Perturbation Theory to cover 
the partially quenched case, explaining how partially quenched Chiral Perturbation 
Theory allows physical results for the Gasser-Leutwyler coefficients to be generated 
from lattice calculations which use the unphysical partially quenched approximation. 

In Chapter ^ our methods for data modeling and statistical error analysis are 
covered. 

In Chapter | we detail the lattice ensembles generated for our study, explaining 
the motivation behind each ensemble's creation. 

In Chapter ^ we present a step-by-step analysis of our lattice data, generating 
values for the Gasser-Leutwyler coefficient combinations 2L^ — L5 and L5 for each 
ensemble studied. A simultaneous analysis of several ensembles, which allows us to 
make a preliminary estimation of the coefficient combinations 2Lq — L4 and L4, is 
also presented. 



In Chapter |10| our final result for 2L^ — L5 is given, along with an analysis of 
our study's systematic error. Using Chiral Perturbation Theory we produce from our 
result a prediction for the light-quark-mass ratio. Secondary results are also given, 
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including values for the coefficients L5, L4, and Lq. The effects of quenching on results 
for the Gasser-Leutwyler coefficients is briefly discussed. 



In Chapter |TT] we summarize our results and discuss the rich potential for future 
work. 
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CHAPTER 2 



QUANTUM CHROMODYNAMICS 



Quarks are the fundamental building blocks of the universe. Bound together by 
gluon exchanges, they are the dominant constituents of hadrons. It is a hadron's 
quarks and the dynamics of the quark-gluon interactions which dictate that hadron's 
characteristics and behavior. 

The physics of quarks and gluons is dominated by the strong force, Quantum 
Chromodynamics (QCD). Working in Euchdean space, QCD is governed by the par- 
tition function: 



where is the gluon field, and q and q are the quark and antiquark fields. The QCD 
Euclidean Lagrangian is: 



The quark field 5 is a vector in three spaces: flavor, color, and spin. The flavor 
index runs from 1 to A'^^, while the color index runs from 1 to iV^ = 3. Finally, quarks 
are Dirac spinors, four component vectors in spin space. 

is the covariant derivative and a matrix in color space, while 7^^ is a Lorentz 
vector of spin space matrices. A4 is the quark mass matrix, a diagonal flavor-space 




(2.1) 



^QCD = ltT[F,,F^''] + q{YD, + M)q 



(2.2) 
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matrix with the quark masses along its diagonal: 



M = diag({m,}) (2.3) 

The gluon field makes its kinetic appearance in the Lagrangian through its 
field strength F^,^: 

= d^A, - d,A^ - ig [A^, A,] (2.4) 

The gluon field is both a Lorentz vector as well as a vector in adjoint color space. 
It has been written as a matrix in color space, with each component of the adjoint 
vector being multiplied by a generator of SU (Nc) : 

A, ^ AIX'^ (2.5) 

The gluon field acts as the parallel transporter of the quarks through color space, 
appearing in the covariant derivative: 

D^ = d^- igA^ (2.6) 

It is through that the quark-gluon interaction arises. 

2.1 Symmetries of Quantum Chromo dynamics 

Of the symmetries respected by the QCD Lagrangian, two of them are of particular 
interest to us: color symmetry and fiavor symmetry, 
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2.1.1 Color Symmetry 

SU (Nc) color symmetry is an exact and local symmetry of QCD under which the 
fields transform as: 

q q' = n{x)q (2.7) 

q q'^qQ^x) (2.8) 

^ A'^ = n{x)(^A^ + ^d,^n^{x) (2.9) 
n(x) = e-'""(^)^" e SU{Nc) (2.10) 

where a°'{x) is a set of A^^ — 1 real functions which parameterize the transformation. 
This transformation of ^4^ results in the gluon field strength transforming as: 

F^,^F'^^^n{x)F,,n^{x) (2.11) 

2.1.2 Flavor Symmetry 

In the case of degenerate quarks: 

M = mql (2.12) 

the Lagrangian respects a global SU{Nf) fiavor symmetry, under which the fields 
transform as: 

q q'^Vq (2.13) 

q q'^qV^ (2.14) 

V ^e-'"""^" e SU{Nf ) (2.15) 

where are the generators of SU{Nf) flavor and v"' are real constants which param- 
eterize the transformation. 



In the case of massless quarks, = 0, known as the chiral hmit, the QCD 
Lagrangian exhibits an even larger SU {Nf)L(^SU (A^/)_r flavor symmetry. The quarks 
split into left- and right-handed pairs: 



q^qL + QR 

Ql = P+q = i(l + 75)9 
qR = P-q^ |(l-75)g 



q^qL + qR (2.19) 

gi = gP_ = ig(l-75) (2.20) 

qR = qP+ = lq{l + j5) (2.21) 
each of which rotates independently under flavor symmetry: 

q q'^LqL + RqR (2.22) 

q ^ q' = qLL^ + qRR^ (2.23) 

L = e-*'"^" e SU{Nf) (2.24) 

R = e"*"""" e SU{Nf) (2.25) 

where and are sets of real constants which parameterize the transformation. 

This symmetry can also be expressed in terms of a vector and an axial vector 
symmetry, SU{Nf)v ® SU{Nf)A- Here the vector symmetry corresponds to the 
smaller flavor symmetry discussed above, where left- and right-handed quarks trans- 
form equivalently, L — R. In contrast, a pure axial vector transformation rotates 
left- and right-handed quarks in opposite directions, L = R^ . The currents associated 
with these symmetries are the non-singlet vector currents: 

j; = q^.T^q (2.26) 



2.16) 
2.17) 
2.18) 



and the non-singlet axial vector currents: 

Jfr = QWr^q (2.27) 

respectively. 

In truth all quarks are neither massless nor degenerate. However, the two lightest 
quark flavors, up and down, have masses and a mass splitting which are quite small 
relative to the typical baryon mass. The next lightest quark flavor, strange, has a 
somewhat small mass, again relative to the typical baryon mass. So, if we restrict 
ourselves to the three lightest flavors of quarks, we would expect SU{Nf = 3)v <8) 
SU {Nf — 3) A flavor symmetry to be a good approximate symmetry of the strong 
interactions. 

However, while SU{Nf)v flavor symmetry is manifest in QCD and its associated 
currents are conserved, SU{Nf)A flavor symmetry does not appear to be respected. 
There exist two possibilities when a theory's Lagrangian contains a symmetry, but the 
theory itself does not appear to respect that symmetry. Either the symmetry is still 
respected but also hidden via spontaneous symmetry breaking, or the symmetry is 
false, ruined by anomalous quantum corrections to the Lagrangian's classical solution. 
The SU (Nf) axial vector flavor symmetry of QCD falls into the flrst category. 

2.1.3 Spontaneous Symmetry Breaking 

Spontaneous symmetry breaking occurs when a theory has not one ground state, 
but rather a set of ground states which transform into one another via the symmetry 
under discussion. In such a situation the theory's vacuum must choose its location 
from among these ground states and thus resides in a position which is not invariant 
under the symmetry. So, while the full theory retains the given symmetry, the state 
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space in the vicinity of the vacuum does not reflect that symmetry. Since this is the 
region important to low-energy interactions and the region explored by perturbation 
theory, the symmetry becomes hidden. 

When the set of vacuum states is continuous, local fluctuations of the vacuum 
within that set are massless. Thus, the theory's spectrum will contain massless parti- 
cles, known as Goldstone bosons, which correspond to those fluctuations. The number 
of spontaneously broken symmetries corresponds to the number of orthogonal direc- 
tions within the set of ground states, and thus corresponds to the resulting number 
of Goldstone bosons. 

It is the interactions of these Goldstone bosons, with one another and with the 
other particles of the theory, which enforces the now hidden symmetry. 

2.1.4 Chiral Condensate 

In massless QCD the energy required to create a quark-antiquark pair from the 
vacuum is small. Because such a pair must have zero total linear and angular momen- 
tum, they will contain a net chiral charge. Thus, the QCD vacuum includes a chiral 
condensate characterized by the non-zero vacuum expectation value of the operator: 



In truth, the vacuum expectation value includes an arbitrary SU {Nf)^ rotation: 



where L — — A e SU{Nf)A- The vacuum is forced to choose the SU{Nf)A 
alignment of this expectation value, thus spontaneously breaking SU {Nf)v®SU (A^/)a 
flavor symmetry down to SU {Nf)v- Note that the chiral condensate is invariant under 



{qq) = {qLqR) + {qRQL) + o 



(2.28) 



(2.29) 
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SU{Nf)y: 



= {(iLV^VqR) + {qRVWqL) = (qq) (2.30) 

where L = R = V E SU{Nf)v- Once the vacuum has made its choice for the 
condensate's ahgnment, we in turn can rotate the definition of our quark fields such 
that the vacuum expectation value does indeed take the form in ( |2.28| ). While the 



freedom for such a redefinition exists in massless QCD, the quark masses of true QCD 
would not allow it. However in that case, the vacuum expectation value will naturally 
align itself with the quark masses, again taking the form in ( |2.28| ). 

Because the set of vacuum states is continuous, connected by elements of SU (A^/)a 
arbitrarily close to identity, we expect to find massless particles in the spectrum, the 
N'j — 1 Goldstone bosons of the spontaneously broken symmetry. The spectrum 
of QCD does not contain any such massless particles. However, together the light 
pseudoscalar mesons — n^, vr^, vr^, K'^, K^, K~^, K~ , and rj — assume the role of the 
Goldstone bosons. Collected, they form an octet of very light particles, one for each 
of the eight generators of the spontaneously broken SU {Nf = 3) a fiavor symmetry. 

It is because the spontaneously broken SU{Nf = 3) a fiavor symmetry is only 
an approximate symmetry of QCD, not an exact symmetry, that these Goldstone 
bosons are not exactly massless. In fact the squared masses of the light mesons are 
proportional to the parameters which break SU{Nf)A fiavor symmetry: the quark 
masses. Thus, the light pseudoscalar mesons are often referred to as pseudo- Goldstone 
bosons. 
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2.2 Uil)A Problem 

In addition to color and flavor symmetry, the QCD Lagrangian is classically invari- 
ant under two additional global U{1) symmetries. The first is vector symmetry, 
under which the quark fields transform by a simple phase rotation: 

q q'^e'^'q (2.31) 

q ^ q' = qe-i^ (2.32) 

This symmetry is exact for arbitrary quark masses and corresponds to the conserved 
singlet vector current: 

J IX = qinQ (2.33) 
d^Jf, = (2.34) 

and to baryon number conservation, a phenomenon clearly observed in experiment. 
The second symmetry is U{1)a axial vector symmetry, under which left- and right- 
handed quarks undergo opposite phase rotations: 

q g' = e*"^5g (2.35) 

q g' = ge^"^^ (2.36) 

Similar to the non-singlet flavor symmetries, this symmetry is only exact in the limit 
of massless quarks. The corresponding singlet axial vector current is: 

4 = QWQ (2-37) 

Assuming degenerate quarks, naive application of the equations of motion results in 
the divergence: 

d^^J^ = 2m,g75? (2.38) 
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which equals zero in the chiral hmit. Thus, we expect the symmetry to be an ap- 
proximate symmetry of QCD. However, U{1)a symmetry calls for degenerate parity 
doublets which are not even approximately manifest in the QCD particle spectrum. 
The mystery of this missing symmetry is known as the U{1)a problem. 

If the symmetry were spontaneously broken, the spectrum would contain a cor- 
responding Goldstone boson. We could imagine combining the singlet axial vector 
symmetry with the non-singlet flavor symmetries: 

U{Nf )A = SU{Nj)a ® U{1)a (2.39) 

and having them spontaneously break together. This spontaneous breaking of U{Nf = 
3) A would result in the eight light mesons mentioned above plus a new ninth light 
pseudoscalar meson. There exists a candidate pseudoscalar meson to fill the roll of 
this ninth light meson: the rj'. However, it has been shown that, if the full U{Nf = 3) a 
symmetry were spontaneously broken, it would constrain the rj' mass: 
01 . The actual rj' is much heavier. U{1)a has no Goldstone boson, and thus sponta- 
neous symmetry breaking can not explain the symmetry's disappearance. We must 
look for a second possibility. 

2.2.1 Adler-Bardeen-Jackiw Anomaly 

When a Lagrangian contains a given symmetry, its classical equations of motion 
will respect that symmetry. However, when the Lagrangian is placed within a path 
integral, that symmetry may be lost. For a given symmetry to survive the transition to 
quantum field theory, not only must the action be invariant under the transformation, 
but the measure of the path integral must also be invariant. When the measure is 
not invariant, the symmetry is said to be anomalously broken. In the context of 
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perturbation theory, this will manifest as a failure of the regularization of radiative 
corrections to respect the symmetry. 

The conserved current corresponding to an unbroken symmetry has a divergence 
of zero. For an anomalously broken symmetry, the divergence of the corresponding 
current will be non-zero. This non-zero divergence is referred to as the anomaly. The 
anomaly which breaks the U{1)a symmetry of QCD is the Adler-Bardeen-Jackiw 
(ABJ) anomaly. 

Under a U{1)a rotation, the measure of the QCD path integral transforms as: 

[VA,][Vq][Vq] ^ [VA,][Vq][Vq]e-''^^^ (2.40) 

where: 

e = ^tr[F^.F'^'^] (2.41) 

and F^iy = \e^,vai3F'^^ ■ This results in the anomalous divergence of the axial vector 
current: 

d^'Jl = 2m,q^5q + O (2.42) 
which no longer equals zero in the massless quark limit. Thus, the U{1)a problem 
appears to be solved. The U{1)a symmetry of QCD is missing because it never existed; 
it is ruined by the ABJ anomaly. However, the situation is not that straightforward. 

2.2.2 Gauge- Variant Axial Vector Current 

The divergence of the axial vector current can also be written as a divergence of 
gauge fields: 

e = ^d^K, (2.43) 



where: 
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= le^uaf^tr [A'd^A^ - '^-igA^A^A^] (2.44) 



We can then define a new current which involves both quark and gauge fields: 




(2.45) 



For massless quarks this new current is then conserved: 



d^'Jl = 



(2.46) 



and once again we find ourselves with a current that appears to be conserved, but 
whose associated symmetry is not manifest in the theory and with which there is no 
associated Goldstone boson. 

The solution to the dilemma centers on the fact that K^j, is not a gauge invariant 



broken, the resulting Goldstone boson might decouple from the theory's physical 
states. When working with a gauge theory in a covariant gauge, the number of 
degrees of freedom in the theory is larger than the number of physical states. A 
condition must be applied to the states of the theory in order to remove unphysical 
states. It is during the application of that condition that the Goldstone bosons of a 
gauge- variant symmetry could decouple completely from the physical states. In fact, 
it has been shown that in QCD they do indeed decouple The modified U{1)a 

symmetry is not observed because it is spontaneously broken, and the Goldstone 
boson which results decouples from QCD's physical states. 

To summarize the fate of the chiral symmetries of QCD, we recall that the massless 
QCD Lagrangian classically respects the symmetry group SU{Nf)v ® SU^Nf)^ ® 
U{l)v ®U{1)a, while only the smaller group of vector symmetries SU{Nf)v ^U{l)v 
is manifest in the theory. The two chiral symmetries have each met a separate fate. 
SU{Nf)A is spontaneously broken by the chiral condensate, resulting in the light 



quantity. Because of this, if the symmetry associated with 



were spontaneously 
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pseudoscalar mesons as Goldstone bosons of this hidden symmetry. U{1)a, on the 
other hand, is gone entirely, forced to gauge variance by the ABJ anomaly. 

2.3 e Vacua 



Generally a quantum field theory which does not experience spontaneous symme- 
try breaking has a single vacuum, the one field configuration with minimum action. 
This is the case when all other field configurations can be smoothly deformed into 
the vacuum state. 

However, if the field theory contains field configurations which can not be smoothly 
deformed into one another, the theory will have multiple vacua. In such a case, the 
space of field states can be grouped into sets of configurations which are smoothly 
connected. In each of these sets, there will be some configuration with minimal action. 
It is these field configurations which are the multiple vacua of the theory. They are 
analogous to multiple local minima in quantum mechanics, separated by infinitely 
high barriers. 

This phenomenon of multiple vacua occurs for certain gauge theories, including 
the SU{N^ = 3) gluon fields of QCD. 

It is possible to write down gluon field configurations which have finite action and 
localized action density, but where the gauge field does not go to zero as we approach 
infinity. This is because a region of space in which the gauge field has the form: 

A^{x) = -n{x)d^n^{x) (2.47) 

contains zero action density. Here Q{x) is an arbitrary smooth function mapping x 
onto the gauge group. As evident from ( p.9|) , the form in ( |2.47|) is simply the gauge 



transform of = 0. When constructing such a field configuration where the action 
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density goes to zero at infinity but the gauge field does not, we must choose fl{x) as 
we approach infinity in each direction. In effect we are choosing a smooth mapping 
of the sphere at infinity Sd-i onto the gauge group, where d is the dimension of the 
Euchdean space. 

As an example, take the case of a U{1) gauge field in d = 2 dimensions. When 
we write down a field configuration, we are choosing a mapping of S*!, a simple circle, 
onto U{1). Since f/(l) is also described by a circle, one is mapping a circle onto a 
circle. Denoting both circles as the phase of a unit vector in the complex plane, we 
can write down examples of such mappings: 

H: e'^ ^ e'^' (2.48) 
The simplest example is the constant mapping: 

Ho : e^'^ ^ 1 (2.49) 
The next simplest is the identity mapping: 

Hi : e''^ ^ e''^ (2.50) 

By visualizing the mappings as vectors located at each point on a circle, we realize 
that it is impossible to deform one of the above mappings into the other using only 
smooth gauge transformations. At some point on the circle, one of the vectors will 
turn clockwise during the deformation, while its neighbor will be required to turn 
counterclockwise. 

A set of mappings which can be smoothly deformed into one another is labeled 
a homotopy class. Thus, the two mappings of ( p.49|) and (|2.5CI|) are representative 



members of two separate homotopy classes. In fact, we can easily write down an 
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infinite number of mappings, each of which belongs to a new homotopy class: 



■V 



e' 



(2.51) 



where v is an integer which enumerates the classes. For any C/(l) gauge field config- 
uration we write down, there is one and only one mapping among those above such 
that the configuration can be smoothly deformed so that its gauge field orientation 
at infinity is described by the mapping. 

The analog for QCD is mappings of SU {N^ — 3) onto 5*3. In this case as well, the 
gauge group is rich enough for there to be an infinite number of homotopy classes. In 
each homotopy class there is a field configuration with the minimum action. These 
are the infinite vacuum states of an SU{Nc = 3) gauge theory. For the homotopy 
class which contains the constant mapping Ho, the vacuum state will correspond to 
the traditional vacuum: a constant zero gauge field. However, for all other homotopy 
classes, the minimum action configurations will contain local kinks in the gauge field. 
These kinks are known as instantons, and each will have a local non-zero action 
density associated with it. 

The winding number u oia field configuration, which is used to label its homotopy 
class, can be shown to equal: 



Thus, we see that the vacuum of QCD is not simple, but rather there is an infinite 
number of potential vacua The true vacuum |^) is then a linear combination of 




(2.52) 
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the In order for this vacuum to have the correct behavior under gauge transfor- 
mations, the coefficients of the combination must have the form: 



|0) = 5^e-V) (2.53) 

u 

where 6' is a new arbitrary parameter of the theory. Incorporating this linear combi- 
nation of vacua into the path integral of QCD effectively adds an additional term to 
the Lagrangian: 

^QCD = 5^e*^^ I [VAMVq][Vq] e-^-^Q^o (2.54) 
= y"[PA^][Pg][Pg] e-^--^QCDcff (2.55) 

where \VA^i, represents a functional integral over field configurations with winding 
number v and: 

-^QCD efl = -^QCD + 

= ^QCD + ^ jl^tr [F^.F'^'^j (2.56) 



From (|2.40|) we can see that the ABJ anomaly has the same form as this new 
term. Thus, if we are assuming massless quarks, we have the freedom to absorb 9 
via a U{1)a rotation of the quark fields. However, for massive quarks, we lose that 
freedom. Because the quark mass eigenstates of QCD are not the same as those of the 
full Standard Model, the actual quark mass matrix M. will include some chiral phase. 
In order to remove this chiral phase and put the mass matrix into the form standard 
for QCD, we apply a U{1)a rotation to the quark fields with a magnitude equal to 
argdet A^. This binds us to that specific rotation, and restricts us from absorbing 9. 
After this rotation, the physical value 9 for the coefficient in becomes: 

^ = ^ + argdetAi (2.57) 
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Each choice of value for 6 corresponds to a unique vacuum choice for QCD. Such 
a continuous set of vacuum states is reminiscent of spontaneous symmetry breaking. 
In fact the 9 vacua are the multiple vacua of the spontaneously broken symmetry 
discussed above. This concurs with the fact that U{1)a rotations move one among 
the ^ vacua. 

2.4 Strong CP 



As demonstrated by ( |2.43| ), is a total divergence and thus can have no ef- 



fect on perturbation-theory calculations. However, it still generates non-perturbative 
symptoms. In particular, breaks CP symmetry and leads to CP violating effects. 
The most significant of these violations is a correction to the zero neutron dipole 
moment. However, the neutron dipole moment is strongly bound by experiment, 
dn < 6.3 X 10^^^ e cm [||. This in turn leads to a bound on ^ 0: 

6 < 10-^° (2.58) 

This extreme smallness of 9 is known as the strong CP problem. 

Were 9 the only parameter in the Standard Model to break CP symmetry, its 
extreme smallness would be of no particular concern. In such a situation, 9 can only 
be multiplicatively renormalized, and thus can remain small over a broad range of 
scales. However, CP is also broken by weak interactions. Thus, even if 9 were small 
at some given scale, at other scales it would be drastically additively renormalized by 
the other CP-breaking elements of the theory. Therefore, the observed smallness of 
9 is deemed unnatural P], and there is likely some as-yet-unknown structure which 
enforces this smallness. As such, the strong CP problem has traditionally been seen 
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as a chink in the armor of the Standard Model, and as a fertile starting point for 
beyond-the-Standard-Model theoretical work. 

One proposed solution to the strong CP problem is through the introduction of 
an additional particle: the axion 0. This field a would couple to the quarks as a 
phase factor on the quark mass matrix: 



^axion = Id^cxd^a + qMe^'^'q (2.59) 



In such a situation both 6 and the chiral phase of Ai can be absorbed via a field 
redefinition of the axion. However, thus far experimental and astrophysical searches 
for the axion have been unsuccessful. 



A second possible solution is the Nelson-Barr mechanism [T^, |TT[. In this sce- 
nario, CP is a symmetry of the fundamental theory and CP violations are due to a 
spontaneous breaking of CP at the GUT scale. Included in the theory are several 
beyond-the-Standard-Model particles, including flavors of heavy fermions. Below the 
GUT scale, 6 gains a non-zero value proportional to the heavy fermion mass divided 
by the GUT scale. The smallness of this ratio, and therefore 6, is no longer unnat- 
ural because, if these fermion flavors were massless, new chiral symmetries would be 
introduced into the theory. In other words, the ratio obtains only a multiplicative 
renormalization as we change scales, and can thus remain small for all scales. 

Another potential solution to the strong CP problem requires that a single flavor 
of quark remain massless. This is in fact the only proposed solution which does not 
require physics from beyond the Standard Model. The obvious candidate for this 
massless quark is the lightest flavor: the up quark. If the up quark were massless, 
we would regain the freedom to apply U{1)a rotations to it, and through the ABJ 
anomaly we could absorb the CP violating term. Surprisingly, current experimental 
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data does not rule out a massless up quark. However, through the lattice calculations 
presented here, we show that a massless up quark is unlikely. 

It should be noted that a massless up quark would not be the end of the dilemma. 
It would simply shift the focus from an explanation of the smallness of 9 to an expla- 
nation of the masslessness of the up quark. However, several beyond-the-Standard- 
Model scenarios for the dynamic generation of quark masses, in which a massless up 
quark is a natural consequence, have been proposed |l^ . 
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CHAPTER 3 



CHIRAL PERTURBATION THEORY 



As a consequence of QCD's non-Abelian nature, it has two important characteris- 
tics: asymptotic freedom and confinement. Asymptotic freedom signifies that at very 
high energies the couphng of the strong force is greatly reduced, and quarks behave 
as if free. In this regime perturbation theory is successful. Conversely, confinement 
indicates that at low energies the interactions are strong enough to confine quarks 
and gluons within bound states. It is for this reason that only color singlet states 
are directly observed, never free quarks or gluons. These low energies are beyond the 
radius of convergence of perturbation theory, and the method becomes useless. 

While it is impossible at low energies to apply perturbation theory to QCD's 
fundamental degrees of freedom, we can use clues from QCD to build an effective 
quantum field theory of its bound states. Such an effective theory is known as Chiral 
Perturbation Theory (ChPT). We focus our attention on to the lightest bound states 
of QCD, the octet of hght pseudoscalar mesons. However, the ideas discussed here 
can also be used to construct an effective theory for the baryons. 

Clues as to the appropriate form for our effective field theory come from QCD 
in the form of symmetries. The Lagrangian of ChPT must respect any symmetries 
respected by QCD. Confinement insures that all bound states are singlets under color. 
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Thus, color symmetry does not restrict the form of our chiral Lagrangian. Both flavor 
and Lorentz symmetry, on the other hand, place strong restrictions on the terms which 
the Lagrangian may include. For the moment we will assume massless quarks, so that 
SU{Nf)v ® SU{Nf)A flavor symmetry is an exact symmetry of QCD. 

Because we have no information other than the symmetries of QCD, we are forced 
to include in our chiral Lagrangian every possible term which respects both flavor 
and Lorentz symmetry. Of course, there are an inflnite number of such terms. So, 
we order the terms by their importance and then ignore those whose importance is 
beyond our chosen sensitivity. The effective Lagrangian is thus an expansion in some 
parameter which establishes the importance of each term: 



9f _ r^(2) , r^(4) , ^(6) 

-Z^ChPT — --2^ChPT + --^ChPT + -^f 



ChPT 



+ 



(3.1) 



Because we are attempting to build a low-energy theory for the light pseudoscalar 
mesons, we will take terms with lower powers of meson momentum to be of greater 
importance. 

We collect the meson flelds tt" into a flavor-space matrix multiplying each 
meson field by its corresponding broken fiavor-symmetry generator: 



$ = ttV^ = (^00 - ^ITr [00] ^ 



1 

71 



du 
su 



ud 

— ^uu + ^dd — |ss 
sd 



us 

ds 



— ^uu — ^dd + 



where: 







u 



d 



= d s] 



(3.2) 



(3.3) 
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<P = TT r = — = 
V2 




71 + 














K- 









The singlet piece of $ has been subtracted, effectively removing rj' from the matrix. 
The fields vr'^ are known as the Cartesian components of the mesons and do not 
correspond to the physical eigenstates of the light-pseudoscalar-meson fields, which 
can be identified as: 

[>° + 7H^ vr+ K+- 

(3.4) 

The unitary matrix S is then built from the meson-field matrix: 

S = e^*'^"""/^ G SU{Nf) (3.5) 
We now define how the meson fields transform under flavor symmetry: 

S ^ S' = LI:R^ (3.6) 
using L and R from (|2.24|) and ([2.251). With this definition $ transforms linearly 



under pure vector transformations: 

$ ^ $' = (3,7) 

where L = R = V. Otherwise, the transformation of the meson fields is non-linear. 
3.1 Leading Order Chiral Perturbation Theory 

Considering all terms allowed by symmetry considerations and then expanding in 
terms of P^/A^ — where p is the meson momentum and is some scale beyond 
which the expansion, and ChPT, breaks down — we find only one meaningful term 
at lowest order: 

-^St = yTrfS.StS^S] (3.8) 
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Note that derivatives of E correspond to powers of meson momentum. While the 
structure of the term is determined by symmetry, its coefficient / is not. The value 
of / is set by the behind-the-scenes dynamics of QCD, which ChPT hides from us. 
In the chiral limit, / equals the pion decay constant: 

/ = A + OK) (3.9) 

where we are using the normalization f^^ ~ 92.4 MeV. Expanding -^chPT terms of 
the meson fields tt" reveals a conventional scalar kinetic term for the mesons as well 
as a tower of interactions involving an increasing number of meson fields. 

Non-zero quark mass breaks the flavor symmetry of QCD, an important building 
block of our chiral Lagrangian. However, there are two critical points which allow 
us to include the effects of quark mass in ChPT. First, we know the form through 
which the quark masses break flavor symmetry: the quark mass matrix Ji4. Thus, 
we can correctly break the symmetry in our effective theory by adding terms to the 
Lagrangian which break SU{Nf)v <8) SU{Nf)A only via insertions of A4. Secondly, 
the quark masses are small, presumably with respect to A^. So, the most important 
of such terms will be those with a low power of M.. 

While the form of Al is known, it is scaled by an unknown constant /i with units 
of mass: 

X = 2nM = 2/idiag({mj}) (3.10) 
In the chiral limit, /j, is directly related to the chiral condensate: 

^ = -7ir + ^W (3-11) 

Now, considering all terms which respect Lorentz and flavor symmetry, except 
via insertions of A4, and expanding simultaneously to 0(p^/A^) and 0(//Al/A^), we 
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determine the LO Euclidean chiral Lagrangian: 



9LO _ (/>{2) 
ChPT — -^ChPT 



(3.12) 




(3.13) 



The low-energy non-perturbative dynamics of QCD have been boiled down, at this 
order, to two unknown constants: / and /i. 

The arrangement of the meson fields in E seems quite arbitrary. In fact there 
are an infinite number of other arrangements and corresponding representations of 
SU{Nf)v <8) SU{Nf)A, each of which would result in a new form for the chiral La- 
grangian. However, any such representation will result in exactly the same meson 
physics as the representation presented here. This is due to a fundamental charac- 
teristic of effective quantum field theories known as universality. For a low-energy 
effective theory, the symmetries of the theory alone determine the resultant physics, 
and the details of the symmetries' representation are unimportant. 

3.2 Leading Order Quark Mass Ratios 

With the LO chiral Lagrangian in hand, predictions can be made for quantities 
in which the constants of the Lagrangian cancel away, such as ratios of semi-leptonic 
decay rates. The meson masses can not be predicted, but their dependence on the 
quark masses can. By inverting those relationships and using experimental values for 
the meson masses, we can calculate the light-quark-mass ratios. 
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The LO ChPT expressions for the meson masses are straightforward: 

Ml, = Ml+ = Ml- = /i (m„ + md) (3.14) 
M|+ = M|_ = /i(m„ + m,) (3.15) 
M|o = M|o = + m,) (3.16) 

These masses Mx are often referred to as the mesons' QCD masses. They are related 
to the physical meson masses by a Quantum Electrodynamics (QED) correction. 
However at lowest order, only the masses of the charged mesons are corrected and all 
by the same amount, a result known as Dashen's theorem [TsU : 

= Ml+-Mlo+0{e^mg) (3.17) 

where Mx represents a meson's physical mass. Thus, QED adds only one additional 
unknown parameter to the picture: 





Mlo 


= /i(m„ + md) 


(3.18) 


Ml+ 


= Ml. 


= /i(m„ + nid) + 6e 


(3.19) 


Ml, 


= Ml- 


= fi{mu + rris) + 6e 


(3.20) 


Mlo 


= M|o 


= fi{md + nis) 


(3.21) 



We are ignoring and will continue to ignore a small correction to M^o due to the 
mixing of the physical 7r° and rj states. The above relationships can be inverted to 
predict the quark-mass ratios [0: 



m„ Ml+ - ML + 2M\ - M\ 



nid Mlo - Ml, + Ml, 

: 20.1 (3.23) 



^ Mlo + Ml, - Ml 
m, Mlo - M2 + M2 
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This LO calculation clearly suggests that the up quark mass is non-zero. However, 
we will see that the NLO calculation muddy the water. 



3.3 Next- To-Leading Order Chiral Perturbation Theory 



When we venture beyond tree-level and attempt to calculate one-loop corrections 
with our LO chiral Lagrangian, we find that the corrections have new matrix struc- 
tures, and thus the divergences can not be absorbed into our LO terms. As it turns 
out, these new terms are the same terms which we would have included in our La- 
grangian had we decided to go out to NLO in meson momentum and quark mass. 
Thus, if we would like to work at one-loop using our LO Lagrangian, we must also 
include tree-level effects from NLO terms. The one-loop corrections will then renor- 
malize the coefficients of our NLO terms. This behavior is a direct consequence of 
ChPT being an effective field theory and will occur order by order. 
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Adding the NLO terms, the Euchdean chiral Lagrangian becomes 



=^ChPT ~ -^ChPT -^ChPT 

= ^Tr [S^St^^S] - ^Tr [St^ + 

- L2 (Tr [S^S+^.S] Tr [S'^S+S^E] ) 

- L3 (^Tr [a^Eta'^Ea.Eta^E] ) 

+ L4 (Tr [a^Eta'^E] Tr [Et^ + xS] ) 
+ L5(Tr[a^Eta'^E(Etx + xS)]) 
-L6(Tr[Etx + xS])' 
-L7(Tr[Etx-xS])' 

- Lg (Tr [EtxEtx] + Tr [x^x^] ) (3.24) 

where Li are additional unknown constants known as the Gasser-Leutwyler (GL) 
coefficients. They are not constrained by chiral symmetry. Rather, they parameterize 
our ignorance concerning the low-energy dynamics of QCD. Terms which couple the 
mesons to a background gauge field have been dropped. Note that each term in 
■^chPT contains either four derivatives, two derivatives and one power of x, or two 
powers of x- 

The magnitude of the chiral scale is an important element in determining 
the validity of ChPT. By looking at the loop corrections to the NLO terms in the 
chiral Lagrangian, we can develop an estimate for A^^,. We find that when the meson 
momentum in loops is cutoff at A = Airf, the radiative corrections from the LO 
Lagrangian are on the same order at the contributions from tree-level NLO diagrams. 
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So, we can only reasonably expect ChPT to be useful below that scale, and thus: 

Ax - 4vr/ (3.25) 

Since / = /tt at lowest order, we find that ~ 1.2 GeV. We are now in a position 
to evaluate the validity of our expansion in the quark mass. At lowest order fiA4 ~ 
M^o, and thus (/i^Vl/A^) ~ 0.2. The quark mass expansion is sound. Therefore, if 
we restrict ourselves to meson momentum below 500 MeV, we can expect ChPT to 
produce accurate predictions. 

3.4 Next-To-Leading Order Quark Mass Ratios 

The NLO expressions for the meson masses are calculable from the NLO La- 
grangian. Using notation similar to the literature, represents the pion mass 
without QED corrections, while represents the kaon mass neglecting the mass 
difference m^^ — m„ and without QED corrections. These uncorrected masses are 
related to the physical meson masses via: 

Ml = Mlo = Mlo (3.26) 
2Mi = M^o + MI+ = Ml, + - (1 + A^;) (M^^ - M^o) (3.27) 

The above equalities are corrected by terms of order O(e^mq) and 0[{md — 171^)^^, as 
well as by NNLO ChPT. The parameter A^; allows for a difference between the QED 
contributions to the pion and kaon masses: 



^Ek 



1 + Ae (3.28) 



^E^ 

Use of Dashen's theorem, as was done in the LO case, is equivalent to using A^; = 0. 



A comprehensive study of experimental data |jT8|, |T9[ suggests a significant deviation 
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from Dashen's theorem and gives the result: 



A^; = 0.84 ± 0.25 (3.29) 
This is the value which will be used in subsequent calculations, unless stated other- 



wise. 



Using the NLO chiral Lagrangian, the uncorrected NLO meson masses are calcu- 



lated In]: 



+ -j^fi{mu + ma) {2L^ - L5) 

16 1 

+ j^l^{mu + md + rus) {2Lq - L4) | 



(3.30) 



M|. = /i(m + m,)|l + 



16 1 

+ —fi[mu + md + ms){2LQ - Li) j (3.31) 



where m = |(mu + m^) and Xx represent chiral logs which arise from divergent loop 
corrections. We choose to cut these loops off at A = A^^ = inf: 



Ml , Ml , , 



The additional correction terms involving GL coefficients come from NLO tree-level 

interactions. 
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We can use these expressions to construct two combinations of the meson masses 



which experience an equal correction Am at NLO |T7| : 



where: 





Mi 


rris + rh 






rrid + rriu 




-Ml, 


_md-mu 




-Ml 


nis — rh 



l + AM + 0(m2) 



1 + Am + Oiml 



(3.33) 
(3.34) 



Am = -X^ +Tn + -pl^i^s - rh) {2Ls - L^) 
= -Z^ +X^ + ^{M],- Ml) {2Ls - L,) 



Thus, a ratio independent of NLO corrections can be defined: 



mz — m 



2 Mi Ml -Ml 
I Ml, - Ml, 



ml - ml Ml Ml,, - Ml 
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Replacing the QCD meson masses with physical masses results in 

Q 



^2 _ 1 Mio + Ml, - (1 + Ae)MI, + (1 + Ae)MI, 



AMI, Mlo - Ml, + (1 + Ae)MI, - (1 + Ae)MI, 

{Mlo + Ml, - (1 + Ae)MI, - (1 - Ae)MIo) 



X 



(3.35) 



(3.36) 



= (22.01 ±0.57)^ (3.37) 

where experimental error in the meson masses is overwhelmed by uncertainty in A^;. 
Use of Dashen's Theorem instead of ( p.29| ) results in a central value of Q = 24.18. 
The quark- mass ratios can be represented as an ellipse [^] which conforms to the 
equation: 

/ m \ I / m \ 

(3.38) 



where m has been taken to be small relative to m„. 
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We now can solve for the quark-mass ratios, but only in terms of the unknown 



NLO correction Am 0: 



2Ml 



m M2(1 + Am) 

u _^ mI mI-{i + /\m)mI 



m. 

771 



(3.39) 
(3.40) 



Setting niu equal to zero and solving ( p. 401) for Am, we find: 

Am = ^ -1 ± Vl + 4Q2 _ 1 

= -0.404 ± 0.002 ± 0.015 ± 0.16 



(3.41) 



where the first uncertainty is due to experimental error in the meson masses, the 
second is due to uncertainty in A^;, and the third comes from an assumption that 
unaccounted for NNLO corrections are on the order of A^. Use of Dashen's theorem 



instead of (13.291) results in the value A 



M 



-0.455 ± 0.002 ± 0.21, for which no 



uncertainly due to A^; is given. The value for Am from ( |3.41| ) corresponds, via 
(Km, to 2^8 - = (-1.25 ± 0.77) x lO^^, or the range: 



-2.02 X 10" 



< 2Ls - L.. < 



0.48 X 10" 



(3.42) 



Using Dashen's theorem corresponds to 2Ls — L5 = (—1.49 ± 0.98) x 10~^. Sharpe 
quotes a range similar to (|3.42| ) in ||22|. The second root of ( p.40|) is not considered 
as it requires the corrections to be even more improbably large. 



From (|3.41|) we see that if the NLO corrections to ChPT are quite large, the 
massless up quark remains a possibility. 



At times, an alternative normalization of the GL coefficients is used: 



(3.43) 
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This is a more natural normalization in that each a, can be expected to be on the 
order of one. Under this normalization the range consistent with a massless up quark 



is: 



-2.6 < 2a8-a5 < -0.6 (3.44) 

As discussed above, one-loop graphs generated by the LO chiral Lagrangian renor- 
malize the GL coefficients. Thus, the GL coefficients are functions of the renormal- 



ization scale. That scale dependence has a simple form |]T7 



L,(A2)=L,(A0 + ^ln^ (3.45) 



where: 



A = - A = — (3.46) 

5 8 ^48 ^ ' 



Li is scale invariant, with = 0. The full set of scaling coefficients can be found 



in [E3i. Unless otherwise stated, we report all GL coefficients at A = 47r/,r 



3.5 Phenomenological Results 

In order to unquestionably rule out a massless up quark, the value of 2Ls—L^ must 
be determined. Because the GL coefficients encapsulate the low energy dynamics 
of QCD, they could in principle be analytically calculated directly from the QCD 
Lagrangian. However, the limitations of perturbation theory at these low energy 
scales has stymied such efforts. Thus, the primary source of information about these 
constants is phenomenological studies. 
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3.5.1 Meson Decay Constants 

Calculating NLO expressions for the meson decay constants using the NLO chiral 
Lagrangian, we find |T^: 



^ = / 1-2J^-J, 



K 



4 



+ -^fx{mu + md + ms)L^~^ (3.47) 



Ik — f \^ — l^TT — 2^K — 4^7? 

4 



+ ^n{mu + + m^) L4 1 (3.48) 



Their ratio isolates 



where: 



1 + Af + Oimi) (3.49) 



4 

Af = - |Jx - iTr, + -pl^ims - m)L5 

= fx. - |Jk - fx, + ^(Mi - M^)L, (3.50) 

We can then fix L5 using experimental values: 

L5 = (0.51 ±0.47) X 10"^ (3.51) 

where the uncertainty comes from both experimental error and an assumption that 

the unknown NNLO corrections to ( p.49|) are on the order of A'j. Use of Dashen's 



theorem instead of ( p.29|) has very little effect on the value, resulting in L5 = (0.50 ± 

0.47) X 10-3. 
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3.5.2 Gell-Mann-Okubo Relation 

At leading order ChPT confirms the traditional Gell-Mann-Okubo relation: 

Ml + ml = AMI (3-52) 

At NLO, CliPT predicts a correction which accounts for the meson masses' observed 
deviation from the relation: 

AMI -Ml -ml , ^ 

AcMO - ' = 0-218 (3-53) 

This NLO correction is: 

AMl,lK-MlI.-3Mllr, 



^GMO 



2- 



Ml - Ml 



- ^ {Ml - Ml) {l2Lr + 6Ls - L,) (3.54) 
Using the observed deviation we can bind a linear combination of the GL coefficients: 



I2L7 + 6L8 - L5 = (-1.08 ± 0.24) X 10-=^ (3.55) 

where uncertainty in and the light meson masses is dwarfed by an assumption that 
the unknown NNLO corrections to the Gell-Mann-Okubo relation are on the order of 
Agmo- Use of Dashen's theorem instead of (|3.29| ) results in the value 12Lj+6Ls—L5 = 



(-1.10 ±0.26) X 10-^ 

With L7 unknown, the combination 2Ls — L5 remains undetermined. 

3.6 Kaplan-Manohar Ambiguity 

We will find that it is impossible to fix the value of 2Ls — using only expres- 
sions from ChPT combined with experimental measurements. The chiral Lagrangian 
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contains an ambiguity among its parameters which prevents us from using the predic- 
tions of ChPT to determine Lq, L7, Lg, and the quark-mass ratios. This ambiguity 
is known as the Kaplan-Manohar (KM) ambiguity pO |. 



The chiral Lagrangian is invariant under a certain redefinition of its coefficients. 
To see this we shift the quark masses in the following way: 

ruu >-> m„ = mu + Xrudms (3.56) 
rrid ^ md = md + Xrriums (3.57) 
nis >— > rhs = riis + XuiuTrid (3.58) 

where A has units of inverse mass. Stated in terms of Xi this redefinition has the 
form: 

X ^ X = X + Ax"^detx (3.59) 

where: 

A = A (3,60) 
Using the Cayley-Hamilton theorem for a 3 x 3 matrix B: 

ldeiB = B^ -B^TiB- \B Tr B^ - |i3 (Tr B) ^ (3.61) 

we can write the shift in x as: 

Xx'^ det X = Ax^Met [x^] 
= a|sxSxS - SxSTr[xS] - iSTrf^S^S] - is(Tr[xS])'| (3.62) 

where det S = 1 has been used. From this we find: 

Tr^x] = Tr^x] - ^jTrfx^xS] - (Tr[xS])'| (3.63) 
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Thus, the mass term in the chiral Lagrangian becomes: 

Tr[Stx + xS] =Tr[Stx + xS] 

+ ^(Tr[Stx + xS])' 

+ ^(Tr[Stx-xS] 
A 



--(Tr[sWx + xSxS]) (3.64) 



From ( ^.241 ) we see that a redefinition of three GL coefficients will absorb these addi- 



tional terms: 





Lq — Lq 


-A 


(3.65) 


Lr >- 


-> Lj = Z77 


-A 


(3.66) 


Ls >- 


Ls = Ls 


+ 2A 


(3.67) 



where: 

A = ^A (3.68) 

Stated explicitly, the KM ambiguity is the invariance in form of the chiral La- 
grangian under the redefinitions ( ^.59| ), ( p3.65| ), ( ^.66| ), and ( |3.67| ). While we had 
believed that there was only one well-defined Lagrangian for ChPT, we now discover 
that there exists a family of Lagrangians, all equally valid and connected via the 
above transformations. 

The KM ambiguity implies that all expressions for physical quantities obtained 
from the chiral Lagrangian will be invariant under the parameter redefinitions. Thus, 
combining these expressions with experimental results can never allow us to dis- 
tinguish between one coefficient set {x^ Lq, L^, Ls) and another (x, Le, I/y, l/s)- For 
example, the combinations of GL coefficients L5 and I2L7 + 6Ls — -L5, which we were 
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able to determine using expressions from ChPT, are both clearly invariant under the 
redefinitions. The ratio defined in ( |3.36| ), which was also fixed using ChPT, is 



invariant to this order under the redefinitions. This becomes clear when we note that 
the squares of the quark masses transform as: 

ml = ml + 2Am„mrfms + Oim^ (3.69) 

Even the meson mass expressions are invariant to the order at which we are working: 



Ml = iiifhu + md) |l + X^ - \lri 

8 

+ j^fJ'irhu + rhd) (21/8 - ^s) 

16 ° 1 

+ j^fi{mu + md + rhs) {2Lq - L4) | 

= Ml + 2Xfimsm + 2fim!^^fim{4:~X) + ^^(2m + m,) (-2A) | + •• ■ 

= M + 2Xfimsm - X—p-fimsm H 

/ 

= Ml + --- (3.70) 



Note that while some measurable quantities appear to break invariance at higher 
order, they in fact continue to be invariant order by order. 

The quark- mass ratios, as well as the quantity 2Ls — L^, are not invariant under 
the redefinition. Thus, using ChPT alone, they can not be fixed. We can only hope 
to determine a one-parameter family of allowed values. 

The KM ambiguity does not represent a symmetry of either QCD or ChPT. It is 
nothing more than an ambiguity in the effective theory's couplings. Thus, true values 
for the quark- mass ratios and the GL coefficients do exist. We have simply found 
that determining those quantities requires theoretical input from outside ChPT. 
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ChPT does place one constraint on Lg, -L7, and Lg. The GL coefficients must 
be of natural order. If they were not, it would imply that ChPT is based on a poor 
expansion, and ChPT would have proven useless. Yet, ChPT's accurate predictions in 
other contexts imply otherwise. We are curious what sort of bounds this naturalness 
places on the up quark mass. If we assume that the maximum reasonable shift in an 
Qij is on the order of one: 



(4:71 fY 

ai = ai + ^-^X = ai + 0{l) (3.71) 
4// 

we find an estimate for the largest reasonable value for A: 

A«-lE^.4^&J,«± (3.72) 

(47r/)2 rus vV^ 

This A corresponds to a shift in m„ on the order of the down quark mass: 

fhu — TUu — Xuidnis ~ nid (3.73) 

Thus, we find that a massless up quark is within the natural range and is not shown 
by this argument to be beyond the scope of possibility. 

In light of the KM ambiguity, we see that we were somewhat naive in our con- 
struction of the chiral Lagrangian. We had stated previously that we knew exactly 
the form of the chiral symmetry breaking structure, and that we could account for 
that breaking via insertions of 7W. Yet, we now find that there exists a continuous 
set of matrices: 

M{X)=M + XM-^detM (3.74) 

each of which is an equally valid choice for breaking the chiral symmetry of our La- 
grangian. If fact, even if the up quark were massless, the form of the chiral symmetry 
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breaking term could naively be mistaken for a non-zero up quark mass: 



MiX) 



mu=0 



ma 



(3.75) 



It is impossible for ChPT to distinguish between the effects of a non-zero up quark 
mass and certain large NLO corrections. 

3.7 Theoretical Estimates 

Various theoretical approximations can be made in an attempt to estimate Am 
and the Gasser-Leutwyler coefficients. 

3.7.1 Resonance Saturation 

While each interaction term in an effective field theory such as ChPT occurs 
at a point, in the full theory they correspond to short-distance interactions, each 
encompassing a tower of graphs involving heavier particles. These particles are heavy 
in the sense that they are more massive than the scale of the effective theory. In the 
case of our ChPT, these heavy particles include all bound states of QCD heavier than 
the light mesons. When the transition is made from the full theory to an effective 
theory, the heavy states are integrated out. This integration shrinks the short-distance 
interaction to a point and condenses the effects of the heavy particle exchanges into 
the coefficient of the effective theory's corresponding couphng term. 

Thus, if we can identifying the most significant of the heavy particle exchanges 
corresponding to a given coupling term, we can estimate the full integration process 
by merely integrating over the identified important exchanges. This leads to a rough 
estimate of the couphng's coefficient in the effective theory. 
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For certain simple coupling terms, the most significant of the contributing interac- 
tions is the exchange of a single heavy particle, the lightest of the heavy particles with 
the correct quantum numbers to mediate the coupling. In such a case the resulting 
estimate for the coupling constant is proportional to the inverse square of the heavy 
particle's mass. 

Both L5 and Lj are examples of coupling constants whose value can be estimated 
as described. We can model their corresponding vertices as being saturated by ex- 
change of the members of the scalar octet and 77' respectively ||2^ , ^ . Adding factors 
which arise from the integration, the estimates become: 

L5 ~ ~ 2.3 X 10-3 Lj ~ ^ = -0.2 X 10-3 (3.76) 

^ AMI 48M3 ^ ' 



with Ms — 980 Me V. Comparing with the determined value of L5 ( |3.51|) , we can see 



that the estimate is off but is of the correct order of magnitude. This allows us to 
approximate the uncertainty of the L7 estimate: 

L7 = (-0.2 ±2.5) X 10-3 (3.77) 



Combined with (|3.51| ) and ( p.55| ), ( p. 77] ) leads to estimates for 2Ls — L5, Am, and 



the light-quark-mass ratio: 



m 



0.4 ±2.1 (3.78) 



We can see that the uncertainly in the resonance saturation estimate is too far sig- 
nificant to rule out a massless up quark. 

3.7.2 Large 

In the large-A^c limit the ABJ anomaly is suppressed and massless QCD has a 
full U{Nf = 3)v ® U{Nf = 3) A flavor symmetry. If we assume that this symmetry 
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spontaneously breaks to U{Nf = 3)v, it would result in nine Goldstone bosons as 



discussed in Section p.2| . A now light r]' would join the ranks of the light mesons and 
take on the role of the ninth Goldstone boson. At NLO in 1/Nc, the anomaly and 
the 7]' mass return. 

The perspective of large Nc can be incorporated into ChPT by constructing a 
Lagrangian which is not only an expansion in P^/A^ and fiAi/A"^, but also in 1/Nc 
| 26| . Additionally, the now light rj' can no longer be excluded from the low-energy 



theory. Instead, we include it as the trace of our meson field matrix $, making S 
an element of U{Nf). This procedure is complicated by the fact that, while the 
ABJ anomaly is suppressed, it is not absent. Correspondingly, the rj' is light, yet 
the symmetry group which must be respected by our Lagrangian is the anomalously- 
broken chiral symmetry, U{Nf)v SU{Nf)A- Under this reduced group, Tr[lnS] oc 
Tr $ is invariant. Thus, the symmetry leaves the Lagrangian's dependence on Tr [in S] 
unconstrained, and the Lagrangian must incorporate arbitrary functions of the rj' field. 
In the end, however, we are saved from these arbitrary functions by the fact that their 
Taylor expansion is an expansion in l/Nc- Thus, large Nc truncates the functions to 
simple forms. 

An analysis of the rj and t]' masses in large- iVc ChPT out to order 0{Ncp'^ /hj^, 
0{Nc^iM/A.\), and 0{N^ = 1) implies a constraint on Am 0: 



Am > ^ = -0.07 3.79 

4Mi-4M2 ^ ' 



Comparing the constraint to (|3.41| ), we can see that laige-Nc considerations suggest 

that a massless up quark is unlikely. 
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3.8 Indirect Phenomenological Results 

To continue our attempt to determine Am via phenomenological results, we are 
forced to draw from evidence beyond the light-meson sector. These results are often 
introduced to light-meson ChPT in the form of the quantity R, the ratio of the 
strength of SU{Nf = 3) breaking over the strength of isospin breaking: 

^ rUs — rh , ^ 

R = — 3.80 

rrid - rriu 



Via (|3.34| ), knowledge of R allows for a determination of Am- 



Results from other systems tend to give somewhat consistent values for R at 
leading order. However, consistency in LO results does not preclude large NLO 
corrections. Additionally, spread in the LO results is significant enough to suggest 
somewhat sizable NLO corrections, perhaps sizable enough to allow for a massless 
up quark. Unfortunately, each of the phenomenological results available requires 
theoretical assumptions in order to tackle their own NLO corrections. Thus, none of 
them provide a model-free determination of the light-meson sector's NLO corrections. 

3.8.1 ip' Branching Ratios 

In the limit of degenerate quarks, the decays: 

^'^J/V^ + 7r° (3.81) 

i)' ^ J/ij + T] (3.82) 

vanish. The first decay is allowed by isospin breaking, while the second is allowed by 
SU {Nf = 3) breaking. Thus, at leading order the ratio of their amplitudes allows for 



a measure of i? [28 



T{^'^J/^ + r]) AR^ ^' ^ ' 
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where Am represents NLO corrections. Attempts have been made to account for the 



NLO corrections resuhing in the value: 

i?^/ = 30 ± 4 (3.84) 
However, the methods are not direct and various theoretical assumptions are required. 



In fact doubts have been raised concerning one of the primary assumptions |30 

3.8.2 p^-u Mixing 

In the limit of perfect isospin symmetry, the uj would be a pure isospin singlet. 
However, isospin is broken and thus the p° and uj mix. The strength of this mixing 
allows for a measurement of R based on the decay: 

uj Tt^-K- (3.85) 



This results in the value pT|: 



i?^ = 41 ± 4 (3.86) 
which does not account for NLO corrections to the vector meson masses. 

3.8.3 Baryon Masses 

At leading order the mass splitting of the baryon octet leads to three independent 
measurements of R. Corrections out to order 0{m?) have been accounted for 



although theoretical assumptions were required. The nucleon, S, and S splittings 
result in the values: 

i?7v = 51 ± 10 (3.87) 

i?s = 43±4 (3.88) 

i?s = 42 ± 6 (3.89) 
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respectively. 

3.8.4 Accepted Values 

Lg has a generally accepted value, often quoted in reviews and first presented in 
|17| . This value is based on the R obtained by averaging the baryon mass splitting 
predictions and the p^-u mixing prediction: 

/2 = 42.6 ±2.5 (3.90) 

where the value used for is more recent than that used in ||T^. This value for R 
results in: 

Am = 0.179 ±0.097 (3.91) 

2L8-L5 = (1.50 ±0.46) X 10-3 (3.92) 

Lj = (-0.55 ± 0.14) X 10-^ Lg = (l.OO ± 0.33) x lO'^ (3.93) 

in 

— = 0.608 ±0.056 (3.94) 



where equations ( p.34| ), ( p.35| ), ( p.40|) , ( p.51| ), and (|3.55|) have been used. The uncer- 



tainty in Am due to NNLO corrections has been assumed to be on the order of A^^. 
Because ( |3.34|) is sensitive to the value of A^;, use of Dashen's theorem instead of 



( p.29| ) results in significantly different values: 

Aj^^ = -0.022 ± 0.057 (3.95) 

2L8 - L5 = (0.55 ± 0.27) X 10-3 (3.96) 

L7 = (-0.31 ±0.12) X 10-3 Lg = (0.53 ±0.31) X 10-3 (3.97) 

ITl 

— = 0.539 ± 0.043 (3.98) 
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Any attempts to calculate rriu/md through a determination of R will encounter this 
sensitivity to the poorly understood QED mass contributions to the light meson 
masses. 

A recent comprehensive study of the relevant experimental data can be found in 
p!9| . The analysis culminates in the result: 

Tn 

— = 0.46 ±0.09 (3.99) 

rrid 

Although the analysis of the data is sophisticated, it still requires the use of indirect 
phenomenological results and various theoretical assumptions in order to determine 
the strength of SU{Nf) breaking. 

All of these results clearly suggest that the up quark is massive. However, the value 
for R was determined using model-dependent assumptions about NLO corrections to 
quantities outside the light-meson sector. We then used this value to calculate the 
NLO correction to the light meson masses. The validity of this process is somewhat 
questionable. 

3.9 First Principles Calculation 

The strong consequences of a massless up quark and the importance of the strong 
CP problem make a first principles calculation of 2Ls — L5, free of model-dependent 
assumptions, very desirable. Currently the only context in which such a calcula- 
tion can be attempted is Lattice Quantum Chromodynamics. Lattice QCD allows 
for a direct and non-perturbative measurement of the Gasser-Leutwyler coefficients, 
numerically evaluating the underlying QCD dynamics from which they obtain their 
values. 
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3.10 Degenerate Quark Masses 

As will be discussed in Chapter ^ in the context of a Lattice QCD calculation, 
we have the freedom to choose the masses of our quarks. For our study, we have 
chosen to use Nf = 3 degenerate light quarks. Thus we present here, assuming 
rriq = m„ = 171^ = and cutting off loops at A = Airf, ChPT's NLO expressions for 
the chiral pseudo-Goldstone boson mass: 

+ ZTJiq (2as — as) + zniqNf {2aQ — 04) j- (3. 100) 



and decay constant: 



A = / i 1 + ^-^lnzmq + zrriq^ + zm^Nf^ \ (3.101) 



where we introduce: 

_ 2/x 



(3.102) 



(4vr/)2 

and we have made use of the alternative normalization of the GL coefficients. Note 
that, other than the fact that we have now left Nf unspecified, these expressions 
follow directly from ( p.3(]| ) and ( |3.47| ). 
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CHAPTER 4 



LATTICE QUANTUM FIELD THEORY 

Lattice Quantum Field Theory (LQFT) is a first-principles non-perturbative nu- 
merical approach to Euclidean-space quantum field theory. 

4.1 Discretization 

It begins with the Euclidean-space partition function of a field theory: 




where the theory contains some set of fields 0^, and -S'f^a] is the Euchdean action for a 
given field configuration [cpa] ■ Any physical observable of the theory can be expressed 
as the expectation value of an operator: the value of the operator evaluated under 
the distribution defined by the partition function: 

{O) = J {UP<l^c?j 0[<P,] e-^[^«l (4.2) 

where the operator 0[4>a] is some mapping of field configurations into real numbers. 

In order to manage the theory numerically, continuous space is replaced by a 
discrete lattice of points, and the infinite extent of space is made finite and compact. 
The number of lattice sites along a given direction /i is denoted by L^, which in all 
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cases we will take to be even, while the distance between lattice sites is denoted by 
a. The theory's fields are attributed values only on the discrete set of locations Xi: 

Ui-f, = a~^Xi-f, G (4.3) 

where each component of is confined to integers in the range [0, L^ — 1]. Derivatives 
within S[(j)a] must also be discretized, ( |A.17| ) and ( |A.18| ). A single field configuration 



now contains a finite number of degrees of freedom, and the functional integral is 
replaced by a finite product of standard integrals: 

(O) =Z 'UU d<P^;n}j 0[<Pa] e-^[<^"l (4.4) 

where Z has been redefined accordingly. In the limit of a small spacing between the 
lattice sites, a — > 0, and a large lattice extent, ^ oo, the results of the lattice 
theory will coincide with the continuum theory. 

While our express motivation for discretization is to allow for a numerical ap- 
proach to quantum field theory, it is worth noting that the discretization procedure 
is also a valid regularization scheme. The finite lattice spacing results in a Lorentz- 
variant ultraviolet momentum cutoff at = 7r/a, removing infinities which arise in 
a perturbative formulation due to loop corrections. Renormalized physical quanti- 
ties become functions of the lattice spacing and remain finite in the continuum limit. 
However, because the regulator breaks Lorentz invariance, it is cumbersome to work 
with and has a narrow range of practical applications. 

While a single field configuration contains a finite number of degrees of freedom, 
the integral in ( [4.4| ) runs over an infinite number of such configurations. Thus, in order 
to attempt the integral numerically, we must apply Monte Carlo techniques. A finite 
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sampling of the infinite field-configuration space is generated, where the probability 
that a given field configuration [0a] is included in the sample is: 



= ^-lg-S[<^a] (45) 

where H^[0a] is the Boltzmann weight of the configuration. Such a set of field con- 
figurations is referred to as an ensemble. In the limit of having a large number of 
configurations in an ensemble, the expectation value of an operator is simply the 
average of its value evaluated on each configuration in the ensemble: 

(^) = ^E(^)W. (4-6) 

n 

where the ensemble contains N configurations and {0)[^^t^^ = C[0a]„ denotes the 
evaluation of the operator on the fixed field configuration [0a] „, the n-th field config- 
uration of the ensemble. 

A given quantum field theory contains some number of fundamental constants. In 
order to fix these constants in the context of LQFT, an equal number of observables 
must be calculated. The results of these calculations are then set to experimentally 
measured values, binding the fundamental constants. From that point on, LQFT 
is predictive. Any additional calculated observables must match experiment. This 
procedure is analogous to choosing the renormalization conditions when implementing 
a continuum regularization scheme. 
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4.2 Ensemble Creation 



Much of the computation time required to numerically implement LQFT is con- 
sumed generating the weighted ensemble of field configurations. Thus, identifying 
efficient algorithm is of paramount importance. 

Any algorithm for generating an ensemble with the proper field-configuration dis- 
tribution will involve the iteration of a two-step update process. This process consists 
of proposing the addition of a configuration to the ensemble, and then accepting or 
rejecting that proposal based on an appropriate probability. Such an algorithm will 
generate the correct ensemble if we insure that the probability of proposing a given 
configuration times the probability of accepting that configuration results in the cor- 
rect probability for that configuration's inclusion in the ensemble. 

The most straightforward update process involves proposing field configurations 
generated randomly, with a distribution that is fiat relative the measure of the par- 
tition function. A configuration [cpa] is then accepted into the ensemble with a prob- 
ability based directly on its Boltzmann weight Vr[0a] = e"'^''^"]. If a random number 
in the range [0, 1] is less than VF[0a], the configuration is added to the ensemble. 

However, the volume of field-configuration space is quite large, and VF[0a] tends 
to be very sharply peaked at a specific set of field configurations. Thus, while this 
procedure generates the correct configuration distribution in a straightforward man- 
ner, it is very inefficient. Much of the computation time will be spent generating 
configurations which are subsequently rejected by the acceptance step. 
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4.2.1 Markov Chains 



We require an algorithm which allows us to primarily propose configurations which 
are in the vicinity of the peak, but which does not skew the probability of each 
configuration's inclusion in the ensemble. 

In order to propose configurations near the peak in H^[0a]) we will generate our 
proposal based on the last configuration accepted into the ensemble. The proposal 
configuration is produced by introducing some change in the previous configuration. 
The magnitude of that change must be small enough that we do not stray too far 
from the peak, but large enough that we can hope to sample a substantial volume of 
configuration space with a reasonably sized ensemble. An ensemble generated using 
such a chain of configurations, each one spawned from the previous, is known as a 
Markov chain. 

Most algorithms for generating Markov chains satisfy two conditions: ergodicity 
and detailed balance. While ergodicity is required for any Markov algorithm, de- 
tailed balance is sufficient, but not necessary, to insure the creation of the correct 
configuration distribution. 

4.2.2 Ergodicity 

As we step along our Markov chain, we must insure that the update algorithm 
does not restrict us to some subset of field- configuration space. Rather, whatever 
process we develop for changing the previous configuration and generating a proposal 
must, given an arbitrary number of iterations, span all of configuration space. This 
condition on the update process is known as ergodicity. 
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4.2.3 Detailed Balance 

To correctly construct our Markov chain, we must determine the proposal-acceptance 
probability which leads to our desired distribution. We begin by noting that, once 
an ensemble has the correct distribution, the configuration densities must be in equi- 
librium. That is, the correct distribution is a fixed point of the update process. 
Therefore, for such an ensemble, the probability of adding field configuration [(pals 
the ensemble, given that configuration [0a] a previous configuration accepted, 

must be equal to the probability of adding field configuration [0a] ^ to the ensemble, 
given that [(pals the previous configuration accepted: 

HIMa - [Mb) = HIMb - IMa) (4-7) 

This condition is known as detailed balance. 

Three distinct factors combine to determine P{[(j)a]A ~^ [^'oIb)'- the probability 
of [(PoIa being the previous configuration accepted into the ensemble Pw{[4>a]A)'j 
the probability of proposing [0a] s, given that [(pa] a the previous configuration 
PpH^PoIa ~^ [^a]^)) the probability of accepting [0a] b, given that [(pa] a the 
previous configuration Pa ([0a] a ^ [Mb]' 

P{[Ua - [<t>a]B) = PwmA)Pp{[MA - [4>a]B)PA{[(t>a]A ^ [Ub) (4-8) 

Equilibrium dictates that: 

Pw{[(t>a]A)Pp{[4>a]A ^ [Ub) Pa^A ^ [Mb) 

- PwmB)Pp{[MB - [<Pa]A)PA{[<Pa]B ^ [Ma) (4-9) 

If we require that the equilibrium point corresponds to an ensemble with our 
desired configuration distribution, then the probability of a configuration being the 
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latest configuration accepted is proportional to its Boltzmann weight: 



PW ilM a) - Z-'W[<l>a]A 



Pw{[cf>a]B) - Z-'W[<f>a]s 



(4.10) 



We have said very little about the process which generates a proposal configuration 
from the previous configuration. However, placing a simple condition on it will allow 
us to determine the acceptance probability. We require that the probability of making 
any given change is equal to the probability of making the reverse change: 



Bound by such a condition, we can easily see that an appropriate acceptance 
probability would be: 



If the proposed configuration has a larger Boltzmann weight than the previous con- 
figuration, it is always accepted into the ensemble. If it has a smaller weight than the 
previous configuration, it is accepted with a probability based on the change in the 
action. 

Such a Markov process, which accepts or rejects a proposal configuration with a 
probability based on its difference in weight from the previous configuration, is known 
as a Metropohs algorithm. 



(4.11) 




(4.12) 
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4.2.4 Autocorrelation Length 



Because each new configuration is generated by making changes to the previously 
accepted configuration, any given configuration in a Markov chain will have similar- 
ities to the configurations which come before it. This correlation between ensembles 
along a Markov chain is known as autocorrelation. 

Due to autocorrelation, a single update step does not sample configuration space 
with the correct weight. Instead, a number of update steps must be performed in order 
to generate only a single independent and correctly weighted sample of configuration 
space. The length that must be moved along a Markov chain in order to go from one 
independent configuration to the next is known as the autocorrelation length. 

We can estimate the autocorrelation length by choosing some observable and 
watching its evolution as we evaluate it on individual configurations along the Markov 
chain. The observable will experience fiuctuations whose frequency can be taken to 
suggest the order of the autocorrelation length. Clearly, the perceived autocorrelation 
length will depend heavily on the observable chosen. Because two configurations are 
only truly decorrelated after an infinite number of update steps, any estimation of 
autocorrelation length will always involve a somewhat arbitrary cutoff decision. 

Reducing the autocorrelation length is desirable, as it allows us to sample a larger 
region of configuration space using a shorter Markov chain. The autocorrelation 
length can be reduced by increasing the magnitude of the changes made when gener- 
ating a proposal configuration. However, such gains can be offset by a reduction in 
the acceptance rate. 
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4.2.5 Acceptance Rate 



In order to minimize the computation time spent generating configurations which 
ultimately go unused, the acceptance rate of a Markov chain process must be kept 
reasonably large. This can be done by reducing the magnitude of the changes which 
are made when generating a proposal configuration, and thus reducing the chances 
of straying significantly from the peak in VF[0a]- However, as discussed above, such 
a reduction causes a corresponding increase in the autocorrelation length. Thus, the 
optimal magnitude of the update step can be elusive, but is generally one which 
results in an acceptance rate of approximately 50%. 

In order to retain a reasonable acceptance rate, but still allow large changes in 
the proposal configuration, update methods which move through configuration space 
along, or nearly along, lines of constant action can be used. Minimizing the change in 
action increases the acceptance rate without ruining detailed balance, and conversely 
allows for larger update steps. By interspersing such constant-action update steps 
with standard update steps, ergodicity is retained. 

4.2.6 Thermalization 

When beginning a Markov chain, one must choose an initial field configuration. 
Because this choice for the head of the chain is arbitrary, and not chosen by the 
Markov process itself, that first configuration, and any configurations which are cor- 
related with it, will not have the correct distribution. Thus, the ensemble will only 
obtain the correct distribution in the limit of an infinitely large ensemble, when the 
effects of the earliest configurations have washed out. 
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This process of decor relation from the initial configuration is known as thermal- 
ization. The number of update steps required by thermalization is obviously closely 
related to an ensemble's autocorrelation length. Since we can only work with finite 
ensembles, configurations generated before the Markov chain is thermalized will be 
dropped from the ensemble. We refer to the point on the Markov chain at which we 
first begin to retain configurations as the thermalization point, with Nt denoting the 
number of configurations dropped. Unfortunately, an estimation of the thermaliza- 
tion point involves the same uncertainties and arbitrariness as an estimation of the 
autocorrelation length. 

4.3 Two-Point Correlation Functions 

A particularly useful observable in LQFT is the two-point correlation function: 

C{x) = {O{x)O{0)) (4.13) 

where the operator 0{x) is some function of the fields local to x and corresponds to 
a set of values for the theory's quantum numbers. The Euclidean-space two-point 
correlation function is analogous to the two-point Green's function of Minkowski 
space. 0{x) creates or annihilates at x every eigenstate of the theory with matching 
quantum numbers, each with an amplitude that is dependent upon the operator's 
exact form. Thus, the correlation function gives the amplitude for creating that 
tower of states at the origin, having them propagate to x, and then annihilating 
them. 

The Euclidean-space states of a quantum field theory are defined to be the eigen- 
states of the time component of the translation operator, also known as the transfer 
matrix. These states decay exponentially with time, each picking up a factor equal 
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to its eigenvalue e~-^* after propagating a distance t in time, where E is the state's 
energy. This corresponds directly to the phase oscillation of a Minkowski-space state 
with a frequency proportional to its energy. Knowing this factor, the tower of states 
created by 0{x) can be made explicit: 

C(f =0,t) = ^^(0|O|n)(n|O|0) e"^"* (4.14) 

n 

where the sum is over all states with appropriate quantum numbers, both single- 
and multi-particle, and (n|C|0) represents the amplitude for the operator 0{x) cre- 
ating the state \n) from the vacuum. The factor (2EnV)~^ comes from a relativistic 
normalization of the states, where V is the spatial volume of a time slice. 

Because 0{x) is local in space, the tower of created states includes states of all 
momenta. We can restrict the states created to those with a specific momentum p by 
Fourier transforming the annihilation operator: 

0{p, t)^J2 e~'^"^0{x, t) (4.15) 

Choosing only states with zero total momentum, we simply sum over all spatial 
positions of a time slice: 

O{p = 0,t) = ^O{x,t) (4.16) 

X 

At zero momentum the energy En of a state is reduced to its mass M„. Thus, using 
a zero-momentum annihilation operator, the correlation function becomes: 



C{t) = C{p = 0,t) = (^O(f,t)O(0,0)^ 

X 

= E^IH(^|0)|V^"* (4.17) 

where the sum is now only over zero-momentum states, and the sum over x has 
canceled with the factor of V~^. Note that restricting the annihilation operator to 
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zero momentum is sufficient. Tfiere is no need to Fourier transform tlie creation 
operator. 

If tfie creation and annifiilation points are well separated in time, and assuming 
that there is an energy gap between the hghtest and second-hghtest states, all states 
will be exponentially damped relative to the lowest energy state 1 1) . Thus: 

limC(t) = -i.|(l|O|0)|V^^* (4.18) 
t-»oo ZMi 

At this point the correlation function has been expressed in terms of only two un- 
knowns: the operator's overlap with the lightest state (l|O|0) and the mass of that 
state Ml. 

Using the techniques of LQFT, we can numerically calculate, for a given set of 
quantum numbers, the two-point correlation function C {t) at large time separations. 
Then, analyzing the t dependence of our result, we can extract the creation amplitude 
and mass of the lowest energy state. By repeating this process for all appropriate 
combinations of quantum numbers, LQFT allows us to non-perturbatively calculate 
from first principles the low-energy spectrum of a quantum field theory. 
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CHAPTER 5 



LATTICE QUANTUM CHROMODYNAMICS 



The study of QCD using the techniques of LQFT is known as Lattice Quantum 
Chromodynamics (LQCD). 



We first describe the discretization of gluon fields. 

5.1.1 Gauge Discretization 

Many of the symmetries of a quantum field theory, the most significant of which 
is Poincare symmetry, are lost during discretization and are only regained in the 
continuum limit. Because local SU{Nc) gauge symmetry is the defining symmetry 
of QCD, the discretization process for the gluon fields will focus on preserving this 
symmetry at non-zero lattice spacing. 

The gluon fields act as a parallel transporter, describing the color transformation 
of the quarks as they move through space. A color vector moving through some path 
C picks up a unitary rotation Uc- 



5.1 Gluon Fields 




(5.1) 



62 



where V denotes path-ordering within the exponentiated integral. In order to retain 
gauge symmetry at non-zero lattice spacing, and to preserve the parallel-transporting 
nature of the gluon fields, we choose to work, not with the gluon fields themselves, 
but instead with unitary color-transformation matrices which we assign to the links 
connecting neighboring lattice sites. These matrices transport a color vector from 
one lattice site to the next, and thus are defined using where the path C is now 

the straight line linking neighboring sites: 



(i-x+afi \ 
ig j dx'A^{x')\ (5.2) 



where U^-n denotes the matrix which transports a vector from the site at x to the site 
at a; + a/i, and jl denotes a unit vector in the direction /i. Correspondingly, f/^.„_^ is 
the matrix which transports a vector from a; to a; — afi. Recall that the four-vector n 
has been defined such that its elements are integers at the lattice sites: 

= a^'^Xf, (5.3) 

At lowest order in a, f/^;„ can be expressed in terms of only the parallel component 
of the gauge field at a point halfway along the relevant link: 



U^-n = exp (^iagA^{x + f /i) H ^ 



(5.4) 



Using link matrices to describe our gauge fields allows for the preservation of 
a discretized local color symmetry. The associated gauge transformation is: 

u^-n u^.^ = ^7„_,.^f/^.„^2|^ (5.5) 

a„ G SUiN,) (5.6) 
where fin is ascribed values only on the lattice sites. 
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5.1.2 Gauge Action 

It is clear that, under the gauge transformation described in (p.5D, any closed loop 
of link matrices will be invariant. The simplest of such loops is a square with sides 
one link in length. This loop in known as the plaquette: 

p — w^^^ = U"^ r/"!" JJ U 
Expressing the plaquette in terms of the gluon field and expanding around its center: 

^ ^y,x ^ ^ t' ^ ^ -y- • • • 

_ ^-aA„{x+^) ^-aAf_i{x+j)~a'^A„Afi{x+^ + ^) 

^ ^aA„{x+^)+a^Ai,Ai,{x+^ + ^)^aAfj,{x+^) _^ . . . 
^ ^-aA^{x+^)-aAf,{x+^)-a^A^A^{x+^ + ^)+Y [A^{x+^),A^L{x+^)] 

^ ^aAv{x+^)+a'^A^Av{x+^ + ^)+aA^{x+^)+Y[-^v{x+^)Ai^{x+^)^ + ■ ■ ■ 
^ ga2(zi,,^^(x+| + |)-Zii,^^(a:+| + |)- _j_ . . . 
^ ^a2(zi^^^(z+| + |)-Z\i,^^(x+| + |)- [^^(x+| + |),.4^(x+| + |)]) _|_ . . . 
_ gia^g-^Mi'Cai+f + |) _|_ . . . 

= 1 + ta^gF^^ix + f/i + f z>) - + f/i + f z>))' + ■ • • (5.8) 

where A^{x) = igA^{x), fl = afi, the matrix identity: 

has been used, and terms higher order in a have been dropped at each step. Thus, 
at lowest order in a, the real trace of the plaquette gives us access to the gauge field 
strength at its center: 

Retr[l-P^,;„] = ^tr[(F^,(x + f/i + fz>))'] +0(0^) (5.10) 
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We are now able to assemble an action which is dependent only on the degrees of 
freedom which remain in our gluon field after discretization, and equals the continuum 
action at lowest order in a. Using 



n 

5^ 5^ 5^ Re tr [l - P^,;„] +0(a2) (5.1i: 



a' 



2 

^2 

n fj, w<fj, 



Finally: 



n /I u<^ 

P = -^ (5.13) 

Because the action is constructed only of closed gauge-link loops, we can be sure that 
gauge symmetry has been preserved. 

5.1.3 Gauge Coupling 



It is evident from (|5.12| ) that, for any lattice calculation which includes only gauge 



fields, there is only one free parameter in the action, the strong coupling constant g. 
We might expect the lattice space a to also appear in the action, but in fact it does 
not. This is because the lattice spacing is not a physical parameter, but rather must 
be thought of as a momentum cutoff for our quantum field theory. Thus, the action 
should depend on our cutoff a only through the renormalization-scale dependence of 
the true parameters of the action. 
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We see then that the strong couphng constant and the lattice spacing are not 
independently free parameters. Once (3 is chosen for a pure-gauge lattice calculation, 
both g and a are fixed, related by the strong coupling's renormalization group equa- 
tion. It is this dependence of the strong coupling on the lattice spacing which allows 
us to conceptualize a continuum limit, and which distinguishes a LQFT calculation 
from a standard statistical mechanics system. 

In more concrete terms, when doing such a lattice calculation, we first choose a 
value for jS. Then we calculate some dimensionful physical observable whose value 
has been determined experimentally. The units in our lattice calculation of the ob- 
servable will appear as powers of a. Thus, by setting our calculated value equal to 
the experimentally measured value, the lattice spacing is fixed. This procedure cor- 
responds directly to the application of a renormalization condition in continuum field 
theory. 



5.1.4 Gauge Updates 

The gauge action (|5.12|) has the very convenient property of locality. That is, any 



given field degree of freedom couples directly only to other degrees of freedom in close 
spatial vicinity. Thus, calculating the shift in action due to some change in the field 
variables requires only information local to that change, and Markov-chain update 
steps can be made with relatively little computational effort. 

In the case of a pure gauge calculation, the update step would involve making 
a change to an individual link in such a manner as to preserve detailed balance. 
Only the plaquettes that contain the changed link need to be recalculated in order to 
determine the shift in the action. The update is then accepted or rejected based on 
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the size of that shift. By doing a number of these individual updates to each hnk in 
the lattice, we can quickly generate a new configuration to add to the Markov chain. 

Unfortunately, the introduction of quarks to our lattice calculation will destroy 
locality, significantly increasing the computation time required for an update. 

5.2 Quark Fields 

We will now bring quarks into our lattice theory. The process will be complicated 
significantly by the fact that they are Dirac spinors, and that the Dirac action is 
linear in momentum. 

5.2.1 Fermion Discretization 

First we define a discretized quark field, attributing it values only on the lattice 
sites: 

5„ = a'^q{x) Qn = a'^q{x) (5.14) 



The powers of a which appear in the definition cancel units contained by the contin- 
uum quark field, making the discretized field suitable for use in numerical calculations. 
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To construct a lattice quark action, we must first discretize the covariant deriva- 
tive: 

Di^q{x) = df,q{x) - igA^{x)q{x) 
1 



2a 
1 

2a 



2 



2a% 



q{x + fi)- q{x - fi) 

+ Ihi^ + A) + - Ihi^ - A) 

1 - igaA^{x + ^)^q{x + fi) 
l+igaA^(x - f)]g(x-/i)| + ••• 

+ ••• 



+ 



(5.15) 



Note that in order to calculate the finite difference at site n, the value of the quark 
field at neighboring sites is first parallel transported to n via the appropriate gauge 
link matrices. In this way the quark-gluon interaction manifests in our lattice theory. 

5.2.2 Fermion Action 



The quark action can now be discretized: 
S^[U,q,q]^ fq{x){YD^ + M)q{x) 



1 f 1 



n ^ 

= EE^"<-[^]5m (5.16) 
n m 

where, if we assume degenerate quarks, the interaction matrix M[U] — M^[U] is: 



(5.17) 



Note that the action includes an implied sum over color, spin, and flavor indices and 
that Mn^rn[U] is a matrix in those three indices, as well as in position. We represent the 
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kinetic term of M„ „j[C/] using D„ „j[C/]. The lattice quark mass thq absorbs one power 
of a and becomes a unitless parameter, suitable for use in a numerical calculation. It 
is related to the dimensionful quark mass by: 

niQ — rhq = arriq (5.18) 

Now including both quarks and gluons in our LQCD action, S{^QQY)[U,q,q] = 
^gW] + Sq[U, Q, q], our partition function is: 

^LQCD = J [VU][Vq][Vq] e-^^l^^ e-^r.m^r.Mr.,m{U]qm (519) 

As a fermion, the quark field is expressed in terms of Grassmann numbers. Because 
there is no simple way to numerically account for the anticommutating property of 
Grassmann numbers, we will integrate out the quark degrees of freedom analytically. 
The result is the determinant of the interaction matrix: 

^LQCD = J [DU] e-^«[^l det M[U] 

= J [VU] e-^i^Q«°[^l (5.20) 

where 

Si^ciCB[U]^Sg[U]-\ndetM[U] 

^ Sg[U]- Tr In M[U] (5.21) 

This determinant, or trace, is over all indices of M[[/]. In the case of M^[U], it is 
over color, spin, flavor, and position indices. Note that if the interaction matrix is 
diagonal in flavor space, as is true for M^[U], the determinant can be factored into 
a product of distinct flavor determinants. 
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5.2.3 Color Symmetry 

The full LQCD action continues to have local color symmetry at non-zero lattice 
spacing: 

Qn (in = ^nQn (5.22) 

Qn q'n = ^n^i (5-23) 

0„ e SU{Nc) (5.25) 

5.2.4 Quark Propagators 

In order to calculate a quark-antiquark correlator, we again must first analytically 
integrate over the fermion degrees of freedom. Via the integration, the quark and 
antiquark insertions generate an inverse of the interaction matrix: 

{qaai;nqbfir,m) = ^LQCD ^ ['DU][Vq][Dq] e-^LQCD[C^W,q] ^^^.^ ^j^^^,^ 

= ^LQCD / m] e-^«[^] det M[U] Mpr^mm (5-26) 

where the correlator's color, spin, flavor, and position indices respectively are shown 
explicitly. We see that M[[/]~^ is acting as the quark propagator. 

Because the interaction matrix is constructed to take quarks into antiquarks, its 
transpose appears as the antiquark propagator: 

= ^LQCD / [m e-^^'""' det M[U] muV)aL,mm (5-27) 
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Inspection of the adjoint of the interaction matrix allows us to express the antiquark 

propagator in terms of the quark propagator: 



M 



N 



m,n+/i ^ n;in^m,n—fi 



e TT r 

n n,m—jx n;n—[i n,m+fi 



(5.28) 



Thus, for this interaction matrix: 



{Qaai;nQbf3j ■,r> 



= ^LQCD / im e-^«[^] det M[U] {j.M'^iU]-^,) 



aain,b/3jm 



(5.29) 



Note that 75 is its own inverse. 

While it is instructive to consider the calculation of a quark-antiquark correlator 
due to its simplicity, it is worth noting that, because it is not a gauge- invariant 
quantity, its expectation value on each configuration includes a random phase and thus 
equals zero after the ensemble average. Quark operators with a non-zero expectation 
value are more complex than the simple two-point quark correlator. Using Wick 
contractions, operators with four or more quark insertions can be reduced to a sum 
of products of quark and antiquark propagators. If the operator is such that these 
products of propagators form closed loops, their expectation value with be gauge 
invariant and potentially non-zero. 
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5.2.5 Quenched Approximation 

The quark fields have introduced the determinant of their interaction matrix into 
the Boltzmann weight. From the perspective of the gauge fields, this determinant 
acts as a non-local interaction, coupling each link to every other link in the lattice. 
While it is still possible to account for this complicated Boltzmann weight in our 
Markov process, the computation time required will be significantly increased. 

In order to avoid this increase, and to return to the simplicity of the pure gauge 
action, an approximation is often made within the partition function known as the 
quenched approximation. The quenched approximation assumes that: 



Thus, the determinant can be accounted for by a simple shift in /3, and the partition 
function reverts to that for pure gauge. Calculation of a quark-antiquark correlator 
under the quenched approximation amounts to the expression: 



where ^qqcd is appropriately defined. 

While the quenched partition function does not correspond to any well-defined uni- 
tary field theory, the resulting system of interactions is often referred to as quenched 



Prom the perspective of perturbation theory, removal of the interaction matrix 
determinant from the QCD partition function corresponds to the removal of quark 
loops from Feynman diagrams. 



det M[U] oc e 



(5.30) 




(5.31) 



QCD (qQCD). 
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5.2.6 Partially Quenched Approximation 

Even when the interaction- matrix determinant is included in the Boltzmann weight 
during a Markov process, a second approximation known as the partially quenched 
approximation is often made. 

The partially quenched approximation arises from the fact that the quark interac- 
tion matrix is a function of the quark mass, M(^^)[C/]. Thus, if we wish to calculate 
the expectation value of an operator at a variety of quark masses, we must generate 
a separate ensemble for each quark mass. However, the creation of a Markov chain is 
computationally expensive, much more so than the calculation of a quark propagator 
under that Markov chain. 

In the partially quenched approximation we generate only a single Markov chain 
using a single quark mass known as the sea, or dynamical, quark mass ms- Observ- 
ables, such as a quark propagator, are then calculated under that ensemble using a 
number of different quark mass values. The quark mass used in the observable is 
known as the valence quark mass mv- Partial quenching is thus an approximation of 
the interaction-matrix determinant at the valence quark mass by that determinant 
evaluated at the dynamical quark mass. 

Calculation of a quark-antiquark correlator under the partially quenched approx- 
imation amounts to the expression: 

{Qaai;nQb/3j ;m) 

= ^p^qcd/ im e-^«[^l detM(^,)[t/] M(^,)[^]„-^\„^,^^.^ (5.32) 
where ^pqQCD is appropriately defined. 
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In a similar fashion to quenching, partial quenching does not result in a uni- 
tary field theory. Nonetheless, the system of interactions which results from partial 
quenching is referred to as partially quenched QCD (pqQCD). 

From the perspective of partial quenching, qQCD is the special case of pqQCD 
in which the dynamical quark mass is taken to infinity. The mass term in the quark 
action Sq[U, q, q] dominates, and the dynamical quarks decouple from the gauge fields, 
affecting the partition function only as an irrelevant constant factor. 

Here we have presented pqQCD as a convenient approximation of full QCD. How- 
ever, we will find that extensions to the concepts behind ChPT will take partial 
quenching beyond a mere approximation, allowing us to calculate physical results 
incalculable within unquenched QCD. 

5.2.7 Fermion Doubling Problem 

In order to determine the particles contained in our discretized quark action, we 
must identify the zeros of the momentum-space action. 
Introducing the momentum-space quark fields: 





(5.33) 



allows for a Fourier transformation of the free-field action: 




(5.34) 
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where: 

^(fc) [1] = ^ ^ 7m sin ak^ + tuq (5.35) 

As in the continuum case, the free-field action is diagonal in momentum space. 
For the remainder of this section, we will assume mq — 0. 

Using an infinitely large lattice volume has resulted in a continuous momentum 
space. However, our discretization of position space has resulted in an action which is 
periodic with respect to momentum. Additionally, momentum space itself is periodic, 
with each component confined to a Brillouin zone of length 27ra~^. We choose the 
range of allowed momenta to be — < — ^72a> such that the maxima of the 
action occur on the limits of the Brillouin zone. 

Each zero of the momentum-space action corresponds to a pole in the propagator, 
and thus a particle of the theory. We expect one particle, the particle corresponding 
to the zero at A; = 0. However, the periodicity of the action results in fifteen additional 
zeros, generating fifteen additional species of quarks. This unexpected proliferation 
of quark species is known as the fermion doubling problem. 

The fermion doubling problem is a direct consequence of the linear nature of the 
Dirac action. Had the action been quadratic, as is the case for scalar fields, the 
additional zeros would have been pushed beyond the Brillouin zone. 

Any of these additional quark species can be shifted to the origin via the redefi- 
nition: 

k^i ^ ~ (5.36) 

(X 

where A is one of sixteen possible binary four-vectors, each of which corresponds to 
a single quark species. We define a binary vector to be a vector whose components 
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take on only the values zero and one. Under this redefinition of momentum: 

^7;fsinafc^ (5.37) 



= i 



where, for each quark species, the Dirac matrices have been redefined in order to 
return M^s [1] to its proper form: 



it = i-)""'^^ (5-38) 

75^ = 7^^3V = 75n(-)''^ (5-39) 

For half of the quark species, 75 has switched sign, resulting in a corresponding switch 
in the definition of chirality. Thus, not only does the fermion doubling problem in- 
troduce unexpected quark species to the theory, but it also hopelessly entangles the 
theory's left- and right-handed degrees of freedom, making impossible their indepen- 
dent rotation. Thus, the definition of a discretized version of chiral symmetry becomes 
impossible. 

The inevitability of fermion doubling for any naive action is demonstrated by the 



Nielsen-Ninomiya theorem and is detailed in ^ Q 
5.3 Staggered Fermions 

There are several ways of handling the fermion doubling problem, each of which 
has its own advantages and disadvantages. The original and perhaps most commonly 



used method is known as Wilson fermions ^6|. A term is added to the fermion action 



which grants the extra quark species a large mass, effectively decoupling them from 
the theory. However, this term also destroys all chiral symmetry. Because of the 
importance of chiral symmetry to our study, such a trampling of the symmetry is 
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unacceptable. Thus, we make use of a second fermion formulation: Kogut-Susskind, 
or staggered, fermions p7|, |38|, |39|], a formulation in which a subset of chiral symmetry 



is retained at non-zero lattice spacing. 

5.3.1 Staggered Action 

The staggered fermion formulation begins with a realization that, through a redef- 
inition of the quark fields, we can diagonalize the fermion action with respect to spin. 
Once each spin component is independent, we drop three of the four components, 
reducing the degrees of freedom within the quark species by a factor of four. Finally, 
the resulting sixteen single-spin-component quark species are collected together to 
form four flavors of Dirac quarks. 

The formulation of staggered fermions makes use of several phase functions, which 
we define now: 

U<IJt U>fJ. 

e„ = n(-)"^ (5.41) 
The staggered fermion fields are defined by a redefinition of the quark fields: 

Xn = r„gn Xn = Qn^i (5.42) 

where r„ is a spin-space matrix, defined by: 

r„^n^M'' = ^"^"^3^^r (5.43) 

Recall that the components of n are integers and that 7^7^ = 1. Thus, r„ traverses 
only sixteen possible values. 
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Substituting x a-nd X iiito the single-flavor fermion action: 



n ^ fj. 



1 



^ ] j 2 ^ ^ ^/u;nXn^n7/i7;ir'n Uj^^.^^Xn+jl ^ |^■,n- fiXn- fi^ "^QXnXr 



1 



where the following identities have been used: 



n— /i 



+ mQXnXn 



(5.44) 



r^r = 1 



(5.45) 
(5.46) 



The Dirac spin matrix 7^ in the fermion action has been replaced by a simple scalar 
phase function 77^;^. Thus, the action is now diagonal with respect to spin, and the 
spin components of the redefined fields have decoupled. We now drop three of the 
four spin components, leaving Xn ^ scalar in spin space and a vector only in color 
space. The result is the staggered fermion action: 



S',[U,X,X] = Y.Y.^'^<rrMX. 



(5.47) 



M: 



(5.48) 



where M^[U] is now a matrix in only color and position. 

5.3.2 Shift Symmetry 

In general, a lattice field theory will have a discrete translation symmetry in which 
the fields are translated multiples of the lattice spacing along directions flush with 
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the lattice. In the continuum hmit this discrete symmetry becomes the continuum's 
continuous translation symmetry. 

In the case of staggered fermions, however, the appearance of the phase factor 
Tj^-n in the action, whose value varies between lattice sites, does not allow for a 
straightforward definition of translation symmetry. Instead, the action has the more 
complex shift symmetry, sometimes referred to as C-shift symmetry, defined by the 
transformation: 

Xn X'n = C,y;nXn+y (5.49) 

Xn x!n^Cu;nXn+0 (5.50) 

where i> points in the direction of the shift. It is easy to verify that this is a symmetry 
of the staggered action through use of the identity: 

r]n;v = C;A (5-51) 

Noting Ci/;nCi^;n+i> = C;i> = 1 revcals that the standard translation symmetries can 
still be defined, yet only for translation lengths which are a multiple of two lattice 
spacings. 

The transfer matrix, which propagates states through time and whose eigenstates 
define the states of a theory, is equivalent to the time component of the translation 
operator. Thus, a staggered transfer matrix which propagates states forward a single 
lattice step can not be defined. Rather, the transfer matrix can only be constructed 
such that it propagates states forward two time steps per application. 

5.3.3 Even-Odd Symmetry 

We define even sites to be those for which summing the components of n results in 
an even number. Conversely, odd sites are defined as those for which the sum is odd. 
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Inspection of M^[U] reveals that, for mg = 0, x on even sites couples only to x on odd 
sites. Similarly, x on odd sites couples only to x on even sites. Thus, the staggered 
fermion action has a. U{l)e <S> U{l)o symmetry, defined by the transformation: 

e^^^Xn even n 

(5.52) 

e*"°Xn odd n 
XnC"'"" even n 

(5.53) 

Xne"^"'^ odd n 

For niQ 7^ 0, only a subset of the symmetry remains, that for which — ag- 
The symmetry can also be expressed in terms ofa?7(l)i®[/(l)g symmetry: 

Xn ^ X'n = e'^"'^'-"'\n (5.54) 
Xn ^ = Xne-'("^+^""^) (5.55) 

where the connection to even-odd symmetry is made via ai = ^{ag + fto) and = 
^{oie — cto)- For niQ ^ 0, the U{1)^ symmetry is broken, and only U{l)i remains. 
This [^(l)e symmetry is the one component of chiral symmetry which survives the 
discretization process and is the reason we use the staggered fermion formulation. 

One other consequence of even-odd symmetry is that the square of the interaction 
matrix M^[U]^ M^[U] couples only even sites to even sites and odd sites to odd sites. 
This becomes useful in simplifying certain calculations. 

5.3.4 Quarks 

The staggered fermion action leaves us with sixteen species of particles, each of 
which has a single spin and flavor component. Collecting these degrees of freedom, 
we define our quark field such that the result is four flavors of Dirac spinors. 
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One option is to make the definition in momentum space [HOl HT], E2|, assigning 



each of the zeros of the momentum-space action a distinct spin and flavor index. 
However, such a procedure leads to a highly non-local definition of the quark field. 
Instead, we chose to define the quark field in position space, constructing the quark 
field Q at a single point by gathering together nearby degrees of freedom in x- 

We divide the lattice into hypercubes. A subset of our full lattice, to which each 
hypercube contributes its lowest member site, makes up a lattice with a spacing of 
2a. We label the sites on this lattice with the index h, such that h/2 G Z^, and 
define our quark field Q such that it takes on values only at the sites of this coarser 
lattice. Each spin- and flavor- component combination of Qh is built from a unique 
linear combination of x the sixteen corners of the hypercube h. The quark field 
{Qh)ai is constructed as a 4 x 4 matrix which mixes spin and flavor space, with one 
spin index a and one flavor index i: 

{Qh)ai = ^ ^{^A)ciXh+A {Qhjia = ^'^Xh+A(X\)ia (5.56) 

A A 

where the summation is over all possible binary four-vectors, and visiting each corner 
of the hypercube h in turn. Note that the normalization of this deflnition varies in 
the literature. 

We will flnd that this spatial separation between the various spin and flavor degrees 
of freedom in our quark fleld breaks a majority of spin and flavor symmetry. Since the 
different spin and flavor degrees of freedom of the quark fleld couple to different gauge 
links, they each experience a distinct gauge environment. Only in the continuum limit, 
where the lattice spacing goes to zero and the spatial separation is removed, will the 
full spin and flavor symmetries of our theory return. 
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5.3.5 Quark Bilinears 



With our quark field defined, we can construct quark bilinears, operators con- 
sisting of local quark-antiquark pairs with a specific spin and flavor structure. It is 
through bihnears that we will access the spin and flavor structure of the staggered 
formalism, using them as the creation and annihilation operators of two-quark bound 
states, our lattice mesons. 

In the continuum we represented bilinears with the notation — qTr'^q, where 
r was some set of spin-space Dirac matrices which determined the bilinear's spin 
structure and was some generator of flavor-space rotations which determined the 
bilinear's flavor structure. 

In the context of staggered quarks, bilinears are represented using the notation: 



where 75 is a matrix which contracts the quark field's spin indices and represents the 
spin structure of the bilinear, and is a matrix which contracts the quark field's 
flavor indices and represents the flavor structure of the bilinear. 

Because the number of quark flavors equals the number of spin components, similar 
bases can be used for the bilinears' spin and flavor matrices. Each matrix takes on 




(5.57) 
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the value of one of the sixteen members of a Chfford algebra: 

' 1 



75 e 



7/. 
7m75 

75 



7/.75 



I 75 



(5.5^ 



The members of a Clifford algebra can be expressed in a manner which parallels 
the definition of r„. Treating S and T as binary four- vectors which encode a bilinear 's 
spin and flavor structure: 



75 



A* 



n 



7. 



7f^7f^7f^7^ 



n 



(5.59) 
(5.60) 



With this notation in place, we can express bilinears in terms of the x and x fields: 

= i 5Z 5Z '^^ [r^75rB7^] Xh+AXh+B (5.61) 

A B 

Thus, a bilinear at a hypercube consists of a linear combination of all possible con- 
tractions of X a^iid X within that hypercube, where the coefficients of the linear com- 
bination depend on the spin and flavor structure of the bilinear. 

In truth (|5.61|) is somewhat oversimplified. The definition of the quark bilinears 
is complicated by the fact that link matrices live in the hypercube over which the 
bilinear is defined. Thus, in order for the quark bilinear definition to have the correct 
color structure, link matrices must be included between non-local contractions of x 
and X- 

Qh {is ® ^ir) '^/i = T Tr [1^751^7^ Xh+Al^A,B;hXh+B (5.62) 



A B 



83 



where UA^B;h represents an equal weighting of all possible shortest gauge link chains 
connecting the corner B to the corner A of the hypercube h. Note that because UA,B;h 
is in general a sum of unitary matrices, it is not itself unitary. 

Returning again to the free-field case, we can solve for the coefficients of the linear 
combination: 

Qh [is ®ir)Qh = ^^^T^'^ \^\ls'^Bl^jr] Xh+AXh+B 

A B 



U+S,B+A-) '-^^^{-r'^^-^Xh+AXh+B 

A B 



= Y,H'''^^'-"Xh+AXh+A+s+:F (5.63) 

A 

where: 

^s,r;A = A - (Cs + Vj^) +3 ■ f]s+r (5.64) 
and we have used the following definitions: 



and identities: 



A-f]B = B-CA (5.67) 



Tv[t\Tj,] =A5a,b (5.68) 

r^r^ = (-)^-^"-+-r^+B r^rt, = (-)^-c-+-rt,^^ (5.69) 

We define the addition of binary vectors modulo two, such that the sum of two binary 
vectors is itself a binary vector. When adding binary and standard vectors together, 
as per the case h + A + S + , the addition of the binary vectors is carried out first. 



From ( |5.63|) we see that, for any given S and combination, only sixteen of the 
sum's coefficients are non-zero. Also, those which are non-zero equal a simple sign 
factor. 
84 



The second term in ^s,r;A can be factored out of the hypercubic sum, and con- 
tributes only an overall sign: 

Qh{ls ® ir)Qh = (-)^-^^+- 5^(-)'^^.-;-x.+AXh+A+5+^ (5.70) 

A 

where: 

¥>'s,^;A = A-{Cs + V^) (5.71) 

For a given bilinear jT^^jr we define the binary four- vector V^ jr = S + J^, and we 
define the distance of that bilinear to be the number of non-zero elements of T>s^yr. 



Inspection of ( |5.63| ) reveals that the non-zero contractions within a bilinear are always 
between corners of the hypercube separated by an offset equal to and thus the 

number of links separating those corners always equals the distance of the bilinear. 

The simplest bilinears are those in which S = They are known as local bilinears 
as they have a distance of zero, and thus all of their contractions are local: 

Qh{ls ® is)Qh = J2i-)'''^^'''''^Xh+AXh+A (5.72) 

A 

In the more specific cases of 75 = = 1 and 75 = = 75, the phase factors 
become quite simple: 



Xh+AXh+A (5.73) 

A 

Qh (75 ® ^5) Qh = Y (^AXh+AXh+A (5.74) 

A 

We can generalize our perspective on the structure (75 ® ^jc) beyond its use in 
the notation of bilinears by viewing it as a position- and color-space matrix which 
operates on x and x fi^ld configurations, also referred to as fermion field vectors. For 
the local case the matrix is diagonal, and simply applies a sign factor to the x a-iid 
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X field at each lattice site. In the non-local case the matrix also swaps the fermion 
degrees of freedom around within each hypercube, shifting the x field at corner A on 
each hypercube to corner A + S + and applying gauge link matrices as appropriate. 
Stated explicitly as a matrix operating on x a^^id x field vectors, (75 ® ^jc) has the 
form: 

{is ® = {-t''^n-t^^'^''''5H+A+S+^,raUA,A^s+:F-,h (5.75) 

where h denotes the hypercube containing n, and A denotes its position within that 
hypercube: 

n 

h = 2 - 

12. 

The first two terms on the right-hand side of (|5.75|) apply a sign factor to x, the second 
of which is position dependent, the third term swaps the x degrees of freedom around 
within each hypercube, and the fourth term applies the appropriate color matrices 
for that movement. From this perspective, a bilinear is seen as the application of 
such a matrix to a fermion field vector which is non-zero only at all corners of a single 
hypercube, and then the contraction of the result with the adjoint of the original 
vector. 

5.3.6 Meson States 

Each bilinear corresponds to a set of values for the spin and fiavor quantum 
numbers of our staggered lattice theory and thus create and annihilate all states 
with like quantum numbers. The lightest state created by a bilinear will be a two- 
quark bound state with appropriate quantum numbers, a staggered lattice meson. 
By determining a correspondence between lattice and continuum quantum numbers, 
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A = n- h 



(5.76) 



we can identify what continuum mesons our lattice mesons become as we move to the 
continuum hmit. 

We choose to label each staggered lattice meson using the name of the hghtest 
Standard Model meson which has corresponding spin quantum numbers. Flavor 
structure takes a backseat to spin structure when labeling our staggered states because 
our staggered theory contains four degenerate light quarks, while the Standard Model 
contains three non-degenerate light quarks. Thus, there is no direct correspondence 
between the flavor quantum numbers of our staggered mesons and the mesons of the 
Standard Model. 

The staggered meson states, each listed with a bilinear having matching quantum 
numbers, include: 



pion: Q{75®^t)Q (5.77) 
rho: Q{li^^T)Q (5.78) 
scalar: Q{l^Cr)Q (5.79) 



where i e {1, 2, 3} chooses the polarization of the rho and the flavor structure of the 
states has been left unspecified. 

The flavor structure of each listed bihnear takes on one of sixteen possible 
values. In the continuum limit the degeneracy of the four staggered quark flavors 
causes these sixteen states to be themselves degenerate. However, at non-zero lattice 
spacing, all but one of the flavor structures becomes non-local, flavor symmetry is 
broken, and the masses of the sixteen states split. 
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For each spin structure, one of the flavor structure choices will correspond to a 
local bilinear. In the case of the staggered pion, that state is: 

local pion: Q (75 ® ^5) Q (5.80) 

At non-zero lattice spacing only this local pion remains light, while the other fifteen 
non-local pions gain a mass contribution dependent upon the lattice spacing. 

Thus far, we have ignored two complications which are inherent in the correspon- 
dence between the bilinear operators and the staggered meson states. 

The first complication relates to the time component of a state's spin. Recall 
that the spin of a particle is defined relative to its momentum four-vector. In the 
case of states with zero spatial momentum, which we create by summing over a time 



slice as described in Section [4.3| , the state's momentum four-vector points along the 
time direction. Thus, the time component of our state's spin is undefined, and there 
exist more bilinears than there are states to create. As a consequence, the bilinears 
Q{is ® ^j^)Q and 0(7475 ® ^J^)Q both couple to the same state. 

The second complication arises because states decay exponentially as they prop- 
agate through time, and thus the use of creation operators, or sources, which are 
distributed in time, such as bilinears, is not straightforward. We will find that each 
bilinear couples to states with two sets of quantum numbers, one set which was naively 
unexpected. 

Looking at the two bilinears (^(75 ^^j^)Q and (5(747575 ® ^4^5^jf)<5 and com- 
paring the sign factors of their x with x contractions, we find the factors the same 
on a given time slice. The difference between their sign factors only arises when we 
compare neighboring time slices. The sign factors of one bilinear will be constant 
between time slices, while the other's factors will flip in sign each time step. 



By defining the binary four-vector /C such that — 7574) we can state the above 
exphcitly: 

V's+K,T+K;A = V's,T;A + M (5.81) 

The only distinction between ^'s+k,t+K;A ^^'^ v's,T;A oscillation along the time 
direction represented by the /I4 term. 

If we choose to use a creation source for our mesons which is limited to a single time 
slice, we will always in effect be choosing a source which is an equal linear combination 
of the bilinears Q{^s ® ^j^)Q and (5(747575 (H) C4:C5^J^)Q- The alternating in time of 
one of the bilinear's sign factors will cause their combined contribution to cancel 
on one of the two time slices of the hypercube. The result of this cancellation is a 
single-time-slice source. 

In an effort to generate only one set of quantum numbers, we might consider 
adding a second time slice to our source. This allows us to use exactly the bi- 
linear Q (75 ® ^jf) Q as our source. However, doing so only suppresses, but does 
not eliminate, the coupling to states having the quantum numbers of the bilinear 
0(747575 ® ^4^5^^^)^- This continued coupling is because states decay exponentially 
as they propagate in time. Thus, in order to overlap only with states of one quan- 
tum number set or another, the exponential decay of states which occurs within our 
source, which is distribute in time, must be accounted for. Yet, such an accounting 
can only be done if we have complete knowledge of the spectrum for the two sets of 
quantum numbers in question. Thus, in practice it is impossible. 

The quantum number sets of the bilinears 0(75 ® Cj')Q Q (747575 'S> ^4^5^:f) 
Q each have an associated lowest energy state, their respective staggered mesons. 
These two mesons are often referred to as parity partners, as they have the same 

89 



quantum numbers up to opposite parity. If, instead of attempting to project out the 
entire tower of states of one of the quantum number sets, we wish merely to project 
out its associated meson state, only knowledge of the time evolution of those two 
lowest-energy states is required. That is, we must know the meson masses ahead of 
time. 

If no attempt is made to account for the time evolution of the meson states, 
and instead the exact two-time-slice bilinear Q{ys <H) Cj=')Q is used as our source, the 
amphtude of the parity-partner meson is suppressed by a factor of tanh ^aM, where 
the two mesons are assumed to have similar masses and M is the average of those 
masses. Note that in the continuum limit, a ^ 0, the suppression is complete. This 
is not surprising as the operators are then local in time. 

For any two parity partners, the positive parity state will propagate through time 
in a straightforward fashion. Its negative parity partner, however, will alternate in 
sign each time step. This unorthodox behavior is a direct consequence of the staggered 
transfer matrix being defined only for time translations of two lattice steps. If we 
consider the states only on every other time step, the sign oscillation of the negative- 
parity state is hidden, and both states behave as expected. 

To summarize the core issue, when the bilinear 0(75®^:f)Q is used to create 
states, it will invariably couple both to states with its own quantum numbers and 
to states with quantum numbers associated with the bilinear (5(747575 ® C4:C5Ct)Q- 
Thus, the direct connection between bilinears and meson states is blurred. 

The local pion is free from this second complication, as it has no parity partner. 
The bilinear Q (75 ® ^5) Q creates only states with the quantum numbers of the local 
pion, and creates no negative-parity states. 
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5.3.7 Quark Action 



Using our definition of quark bilinears, we can express tlie free-field staggered 
fermion action in terms of Q and Q: 



1 



Xn+fl Xn- 



n ^ fj, 



h A 



AXh+A 



Xh+A+fi — Xh+A-fi 



+ rriQXnXn 

+ mgXh+AXh+A 



E { i E [^^^ (7,. ® 1) {Qh+2f. - Q h-2fi J 



+ Qh{i5 ® ^M^s) - 2Qh + Qh-2fi)\ + mQQh{i ® i)Q, 

wliere we liave used tfie following identities: 

Xh+A = Tr [QhTa] 

Xh+A+fi = Tr {So^A^Qh + Si^A^Qh+2il)TA+f, 

Xh+A-fi = Tr (5o,A^Q/i-2A + 5i,A^Qh)T\_f^ 

J2 Tr [QhTA] Tr [g.r^r^^r^] = 4Tr [Q^r^g^r^ 

A 

(^A^A = 75^Al5 5o,A^ = 1(1 + (^Ar]f^;ACf^;A) 

Vl^;ACt^;A^A = 1^^^A'1,, = - £^^7/.;^^;^) 



(5.82) 



(5.83) 
(5.84) 
(5.85) 

(5.86) 
(5.87) 
(5.88) 



Tlie first term in 



is expected. It is a discretization of the standard flavor- 



diagonal kinetic term of the continuum quark action. However, the second term 
corresponds to the lattice spacing a times a discretization of the second derivative of 
the quark field ( |A.18|) . Thus, this second term is a lattice artifact. In the continuum 
limit, a — i> 0, it disappears and the staggered quark action becomes the standard 
continuum action. 
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Yet at non-zero lattice spacing, this term remains, and as the term is not diagonal 
in flavor space, it demonstrates explicitly the flavor symmetry breaking present within 
the staggered formalism. 

5.3.8 Flavor Symmetry 

A continuum four-flavor massless QCD theory classically has a U{4)y U{A)a 
flavor symmetry. However, in our staggered discretization of that four-flavor theory, 
the second kinetic term in the staggered quark action ( |5.82| ) breaks much of the flavor 



symmetry, and at flnite lattice spacing only a remnant is preserved. 

For massive quarks, the staggered quark action has a U{l)i symmetry, deflned by 
the transformation: 

Q ^ g' = e-i(i«i)g (5.89) 
Q _ Q' = Qe-*«i(i®i) (5.90) 

For massless quarks, the symmetry expands to U{l)i^U{l)^^, where U{1)^^ is deflned 
by the transformation: 

Q Q' = e'^'^^^'^^'^Q (5.91) 

Q g' = ge*"^(^«®«^) (5.92) 

This U (1)1 ^7(1)75 symmetry is equivalent to the ?7(l)i Cg) f/ (l)e symmetry described 
by ( ^.54|) and (|5.55|) , simply reexpressed in terms of the staggered quark flelds. In 



the continuum limit, the U{l)i symmetry becomes the flavor-singlet vector symmetry 
U{l)v- The U{l)e symmetry becomes a single U{1) subgroup of the continuum's flavor 
non-singlet axial vector symmetry. Following its role as a non-singlet axial vector 
symmetry, f/(l)e is spontaneously broken, with the resulting Goldstone boson being 
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a two-particle bound state created and annihilated by the bilinear 0(75 0^5) Q. When 
we move away from the chiral limit, U{1)^ becomes an approximate symmetry, and 
the bound state becomes a pseudo-Goldstone boson with a squared mass proportional 
to the quark mass. 

The form of the creation bilinear reveals that the Goldstone boson of the spon- 
taneously broken U{1)^ symmetry is the local staggered pion. It is clear now that 
the local pion remains light at non-zero lattice spacing because a remnant of the 
continuum's non-singlet axial symmetry is preserved. It is also not surprising that 
the fifteen other non-local pions gain a mass, as the staggered formulation breaks the 
flavor symmetries for which they would have been Goldstone bosons. 

It is the existence of this robust Goldstone boson, still present an non-zero lattice 
spacing, which will prove to make the staggered fermion formulation a valuable tool in 
our study. Because the local pion is a Goldstone boson, Chiral Perturbation Theory 
allows us to calculate its mass and decay constant in terms of certain GL coefficients. 
By evaluating these quantities using a lattice calculation, we can determine the value 
of those GL coefficients. 



5.3.9 Nfj^4: 

One remaining stumbling block is that the number of quark flavors is restricted 
by the staggered fermion formulation to a multiple of four. Yet, the Standard Model 
contains three light flavors. We resolve this issue by taking the interaction-matrix 
determinant to a fractional power. 

As discussed in Section |5.3.7| , the staggered interaction matrix M^[U] becomes 



flavor-space diagonal in the continuum limit. Thus in that limit, its determinant can 
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be factored into the product of four equivalent determinants, one associated with each 
quark flavor: 



If we desire some number of quark flavors Nf other than four, we need only to include 
Nf powers of the fourth root of the interaction-matrix determinant in our partition 
function. In terms of the effective gauge action, this choice of flavors becomes a simple 
factor: 



Of course at finite lattice spacing, the staggered action is not flavor diagonal and 
its determinant will not factor. Thus, this procedure for choosing Nf clearly becomes 
invalid. Yet, the inclusion of the Nf/A factor in a numerical calculation proves to be 
straightforward and, whatever effect that factor has at non-zero lattice spacing, as 
a becomes small that effect will become an improving approximation of Nf flavors. 
Thus, in our lattice calculations we will use this method to approximate Nf = 3 quark 
flavors, the same number of hght flavors as is present in the Standard Model. 

5.4 Conjugate Gradient 

At their core most lattice calculations involve the inverse of the fermion interaction 
matrix. In the case of our study, the inverse is required during both the generation 
of our ensemble configurations and during the calculation of our meson propagators. 

The staggered fermion interaction matrix M^[U] is very large, L^N^ x L'^N^, but 
since it only connects each site to its neighbors, it is also quite sparse, with only 




(5.93) 



aS 



'LQCD 




(5.94) 
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L^Nc{4:Nc + 1) non-zero elements. In both cases, we use to represent L1L2L3L4. 
Because of this sparseness, the memory required to store the interaction matrix grows 
only linearly with the lattice volume. Its inverse M'^[U]^^, however, is not sparse. 
Thus, we make no attempt to calculate and store the inverse interaction matrix 
directly. Rather, we use a numerical algorithm which allows us to apply the inverse 
matrix to a vector via repeated application of the original matrix. This algorithm, 
which finds use in a broad range of fields, is known as the Conjugate Gradient (CG) 
method [Q. For a clear yet complete explanation of CG, we refer the reader to 



We use CG to apply the inverse interaction matrix to some field vector W, calcu- 
lating X, where X has the form: 

X = M^[U]-^W (5.95) 

In practice however, we can not use CG to apply M^[U]^'^ to a vector directly be- 
cause M^[U] is not positive definite, a characteristic required by CG. Yet the square 
of the interaction matrix M"^[f/]^M'^[f/] is positive definite. Thus, we are able to 
calculate X by first applying M^[U]'' to W, and then using CG to apply the inverse 
of M^[U]^M^[U]: 

X = {M^[U]'^M^[U]y^M^[U]'^W (5.96) 



Recalling Section ^.3.3| , the calculation is simplified by the fact that the square of 
the interaction matrix, and thus its inverse, couples only even sites to even sites and 
odd sites to odd sites. The even components of X can be calculated by inverting a 
matrix which consists of only the even-to-even elements of M^[U]^ M^[U], a smaller 
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matrix which is thus easier to invert: 



Xe = {m''[u]^m'[u])-^'^{m'[u]^w)^ 

= {M'[U]^M'[U]y^^^{-D'[Ul^^Wo + mQWe) (5.97) 

where the subscript e (o) denotes the even (odd) half of the lattice sites. We have 
used D[U]^ = -D[U]. From (|5]9|) we know that: 



Wo={M^[U]X)^ 

= D'[Ul^^X, + mQX, (5.98) 

Thus, the odd elements of X can be put in terms of its even elements: 

X, = ^{Wo-D'[Ul^^X,) (5.99) 

We see that both the even and odd components of our resulting field vector X can be 
computed by inverting only the even half of the interaction matrix. This procedure is 
known as preconditioning and will significantly speed up our numerical calculations. 

5.5 R Algorithm 

The advent of fermion fields in our lattice field theory results in a non-local effec- 
tive gauge action: 



SIqcbIU] = S,[U] + ^TT\nM^[U] (5.101) 

where the Tr In M'^[[/] term mixes all gauge degrees of freedom. This non-locality sig- 
nificantly increases the computational effort required to generate an ensemble. When 
producing a Markov chain using the local pure-gauge action, determining the change 
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in action due to some local modification requires only information in the neighbor- 
hood of the modification. With our non-local action, however, determining the change 
in action due to any modification requires a full recalculation of the configuration's 
action. As such, we can no longer afford to grow our Markov chain via a large number 
of small localized changes. 

Instead, we require a procedure which permits large steps through configuration 
space during which all gauge field degrees of freedom are updated simultaneously. 
These update steps must change the equilibrium ensemble distribution either not at 
all or very slightly, such that there is either no chance, or only a very small chance, 
of rejection. Additionally, the algorithm must permit simulation of Nf = 3 flavors of 
dynamical fermions. 



The R algorithm also known as Hybrid Molecular Dynamics (HMD), fits all 
of these criteria. The strong point of HMD is that it chooses the next configura- 
tion in our Markov chain with the correct probability distribution. No acceptance 
step is required. Additionally, the simulation of any number of dynamical flavors is 
straightforward. Unfortunately, the effective probability distribution used by the R 
algorithm is only accurate to within a certain uncorrectable error. 



5.5.1 Hybrid Molecular Dynamics 

The first conceptual step towards HMD is the introduction of an auxiliary field 
H^-n to the theory. This field consists of traceless Hermitian color matrices, four of 
which live on each lattice site. They are parameterized as: 



a 



(5.102) 
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where A" are the generators of SU (3) color. Being an auxihary field, 11^.^ has a trivial 
action: 



■s _ 

HMD — 




(5.103) 



(5.104) 



(5.105) 



As there are no interaction terms between the auxiliary field H^-n and the gauge field 
Uf^-n, its introduction into the theory does not affect the expectation value of gauge or 
fermion field operators. The quark and gluon physics of the HMD partition function 
is the same as that of the LQCD partition function. 

With its simple action, updating the auxiliary field is effortless. A heat-bath up- 
date can be used, in which the next field configuration is simply chosen using 
the proper probability distribution. This field configuration choice is made by set- 
ting each to a Gaussian random number with appropriate normalization. No 
acceptance step is required. 

Given that [C/]^ was the last gauge configuration added to our Markov chain, after 
including a heat-bath updated auxiliary field [i?]^, we have the HMD configuration 
[H, U]^. What we now require is a procedure which will take us from [H, U]^ to a 
new configuration [H^U]^ with equal HMD action S^y^^lH, [/] in a deterministic and 
reversible manor. Once [H, U]^ is found, [U]^ can be added to our ensemble with no 
acceptance step, as there has been no change in the HMD action. The procedure used 
by HMD to generate [H, U]^ from [H, ^7]^ is based on classical molecular dynamics. 
This is from where HMD draws its name. 
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In the context of HMD, the gauge field configuration [[/] is thought of as denoting 
the position of a classical particle in the very-large-dimensional configuration space. 
The LQCD action 5'fQ(-;j-)[f/] is taken to be a static potential through which the particle 
moves. Finally, the auxiliary field configuration [H] acts as the particle's canonical 
momentum. With these identifications, the HMD action 5'j^MD[-^)^] same 
form as the one-particle system's Hamiltonian M'. 

The configuration [i/, f/]^ is treated as the initial condition for a trajectory through 
configuration space. Using the Hamiltonian's classical equations of motion, we move 
the configuration along the trajectory for some interval in artificial HMD time. By 
the time we reach the endpoint of the trajectory [if, ?7]^, the gauge field degrees of 
freedom have changed significantly. We have taken a large directed step through con- 
figuration space. Additionally, [if, U]^ has the same HMD action as our starting point 
[ii, f/]^. This is because, along a classical trajectory through a time-independent po- 
tential, a system's Hamiltonian is conserved. Thus, we can safely add the gauge 
configuration [f/]^ to our ensemble with no acceptance step. In a sense, when a 
heat-bath update is applied to the auxiliary field, HMD is choosing, with the proper 
distribution, our ensemble's next gauge configuration. It is then simply a matter of 
calculating, via the classical equations of motion, which configuration the R algorithm 
has chosen. 



5.5.2 Ensemble Equilibrium 

To see clearly that this sort of update process is valid, we can investigate the en- 



semble equilibrium condition discussed in Section |4. 2. 3| . That is, an ensemble with the 
desired configuration distribution must sit in equilibrium during the update process: 
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PiiUU - [U]b) = - [U]a) (5.106) 

As before, P{[U]j^ [U]b) has three factors: 

PiPU - [U]b) = Pw{[UU)Pp{[UU - [U]b)Pa{[UU - [U]b) (5.107) 

where: 

PwiiUU) = ^L(icDe-^''^«°'''l- (5.108) 
and as we have asserted that no acceptance step is required: 

Pa{[UU^[U]s)^1 (5.109) 

The remaining factor Pp{[U]^ — > [17]^) is essentially the probability that the heat- 
bath update of the auxiliary field will choose a certain trajectory, the trajectory 
which takes us from the position [U]^ in configuration space to the position [U]^. 
Conversely, Pp{[U]^ [U]a) is the probability of choosing that same trajectory in 
the reverse direction. Associated with the trajectory is a value for the conserved 
Hamiltonian Jf which is independent of the direction the trajectory is traversed. 

The probability of choosing a given trajectory is based solely on the kinetic energy 
Sh[H] required to take the trajectory. The kinetic energy must make up the difference 
between the trajectory's Hamiltonian M' and the static potential at the trajectory's 
starting point SI(^q^[U]^: 

SH[H] = Ji^-SlQCD[U]A (5.110) 

Thus, the probability of choosing to take the trajectory in the forward direction over 
choosing to take that same trajectory in the reverse direction is: 



Pp{[U]b ^ [U]a) e-^+^^QCD[^l« 
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e^lQCDPU-sEacDMB (5.111) 



This is exactly the probabihty required to generate the desired configuration distri- 
bution: 



^Pw{[U]B)PpmB^[U]A) 

= pmB - [u]a) 



[U]a) 

(5.112) 



5.5.3 Finite Step-Size Error 



We numerically evolve an auxiliary and gauge field configuration along its trajec- 
tory using finite steps in HMD time of some chosen size. Preceding each step, the 
gradient of the static potential -S'fqQofC/] is calculated in order to determine the mag- 
nitude and direction of the step. Each such calculation requires an application of the 
inverse interaction matrix, the most numerically demanding aspect of the algorithm. 
As such, it may be tempting to use a large time step in the interests of efficiency. 
However, the finite size of these steps introduces error into the R algorithm's ensem- 
ble distribution. The finite steps cause the numerical evolution to stray from the true 
trajectory. Thus, the end-point configuration will not be exactly the configuration 
chosen by the heat-bath update of the auxihary field. We must insure that the step 
size used is sufficiently small such that our ensemble's configuration distribution is 
not significantly skewed. 
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5.5.4 Nfj^4 



When simulating Nf = 4 dynamical staggered quarks, the $ algorithm [55|, also 
known as Hybrid Monte Carlo (HMC), can be used. This algorithm uses molecular 
dynamics in the same fashion as the R algorithm. The difference between the two lies 
in how they handle the interaction-matrix determinant in the effective gauge action 

'S'lqcd [U] ■ 

In the case of the $ algorithm, the determinant is replaced by a set of pseudo- 
fermion fields whose interaction matrix is the inverse of the staggered quark interac- 
tion matrix. Such a pseudo-fermion interaction matrix generates the appropriate de- 
terminant factor in the ensemble average. As the pseudo-fermions are not Grassmann 
fields, we can work with them numerically, and as their action is straightforward, they 
are easily updated at the start of each trajectory using a heat-bath method similar 
to that used for the auxiliary field. 

A particularly noteworthy feature of the $ algorithm is that the use of the pseudo- 
fermion fields allows the HMC action to be known exactly. As a result, the error 
introduced by finite step size can be accounted for with an acceptance step at the 
end of each trajectory. While the numerical evolution will stray from the correct 
trajectory before it reaches its endpoint, we can make up for this error by accepting 
or rejecting the entire trajectory based on the difference between the the HMC action 
at the beginning and end points of the trajectory 

In the case of the R algorithm, because we are taking the fermion interaction- 
matrix determinant to a fractional power in order to simulate at Nf ^ 4, the deter- 
minant can not be replaced with a functional integral over pseudo-fermionic degrees 
of freedom. Instead, the determinant is evaluated using a noisy estimator. Such a 
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procedure does not allow for an exact determination of the HMD action, and thus 
the error introduced by finite step size can not be accounted for via an acceptance 
step. 

In summary, the R algorithm allows us to take large directed steps through con- 
figuration space with no chance of losing numerical effort due to an acceptance test. 
Additionally, it allows us to operate at Nf ^ 4. Its primary disadvantage is that the 
algorithm is inexact, as the finite step-size errors can not be accounted for. 

For the specifics of the R algorithm, including the exact form of the classical 



equations of motion, the reader is referred to |45 . 



5.6 Hypercubic Blocking 

In order to gain access to the GL coefficients, we are studying the local pion 
mass's dependence on the quark mass. At NLO in ChPT, the local pion propagator 
includes graphs with non-local pion loops. Thus, its mass gains a dependence on the 
masses of the non-local pions, and through them a NLO dependence on the strength 
of the staggered fermion formulation's flavor symmetry breaking. Because this is the 
same order at which the GL coefficients appear, flavor symmetry breaking has the 
potential to introduce significant systematic error in the values we observe for the GL 
coefficients. 

If the lattice spacing of our ensembles were small enough, flavor symmetry break- 
ing would become negligible. However, as is always the case in lattice calculations, 
due to limited computer resources our lattice spacing is not as small as we would like. 
Thus, some other method for controlling flavor symmetry breaking is required. 
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An oft-used method for reducing lattice artifacts is the smearing of gauge hnks. 
Each of a configuration's gauge hnks is replaced by a linear combination of itself and 
other local gauge paths which connect the same endpoints. The result is a smeared, 
or fat, link configuration. Such a process does not affect long distance behavior, but 
reduces lattice artifacts by smoothing out short- distance gauge fluctuations. A reduc- 
tion in statistical error, especially in quantities sensitive to short- distance fluctuations, 
is often seen, as well as an improvement in rotational symmetry. 

In the context of staggered fermions, flavor symmetry breaking is the result of the 
various corners of a hypercube experiencing distinct gauge backgrounds. Smearing 
the gauge links reduces these differences, and thus reduces flavor symmetry breaking. 
This has been demonstrated in both quenched [^] and full LQCD |^, and is also 



supported by perturbation theory |Q. It is not surprising that a procedure which 



in general improves rotational symmetry would, in the context of staggered fermions, 
improve flavor symmetry as well, as the Lorentz and flavor symmetries have been 
tightly entangled. 

Often multiple iterations, or levels, of a smearing process are applied to an en- 
semble, as the effects of smearing become more pronounced with the number of levels 
applied. However, a large number of smearing levels results in a highly non-local 
definition of the fat link and significantly alters short distance behavior. Thus, there 
is a trade-off between using a large number of smearing levels to greatly reduce short- 
distance gauge noise, and using a moderate number of smearing levels to preserve 
short- distance physics. 

From the perspective of reducing flavor symmetry breaking in staggered fermions, 
the optimal point in this trade-off is a smearing procedure which maximally smears 
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a gauge link without reaching beyond that hnk's hypercubes. A procedure developed 



for this exact purpose is hypercubic blocking |50[. It constructs a fat link by 
smearing it only with links contained by the eight hypercubes of which that link is 
an edge. 

In the application of a standard smearing level, a link is replaced by a projected 
linear combination of itself and the six three-link gauge paths, known as staples, 
which connect the same endpoints. With hypercubic blocking, three levels similar to 
this manner of smearing are applied. The defining feature is that staples which reach 
beyond the hypercube of the resulting fat link are not used. As such, the staples 
used in a given smearing level must be orthogonal to the fat link and to the staples 
used in all previous levels. With only three directions orthogonal to the original link, 
hypercubic blocking stops after three levels. 

Taking f/^.„ to be the original thin links of a configuration and V"^;„ to be the 
resulting fat links of the hypercubic-blocked configuration, we layout the definition of 
hypercubic blocking explicitly: 



proj 

SU{3) 



5(7(3) 



SU{3) 



Pf^p. ± 
pi^v 



(5.113) 
(5.114) 

(5.115) 



where = Xl;.„_-. The parameters ai, 0:2, and 03 are the adjustable parameters 

of the procedure. Throughout our work we use values for these parameters suggested 
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by early hypercubic blocking literature 



ai = 0.75 02 = 0.6 03 = 0.3 (5.116) 

Because the sum of a set of unitary matrices is not itself unitary, the result of 
each level of smearing must be projected back into the space of unitary matrices. 
This projection process must result in a unitary matrix which is as close as possible 
to the original non-unitary matrix so that the resulting fat link characterizes as well 
as possible the sum over local paths. It is not obvious, however, what metric this 
measure of closeness should use. As such, we simply define a reasonable method for 
identifying the closest unitary matrix, and do not trouble ourselves with the exact 
form of the resulting metric. For the definition of the SU{3) projection and details 
of the procedure used to implement that definition, see Appendix 0. 

It has been demonstrated that hypercubic blocking reduces the mass splitting 



between the local and lightest non-local staggered pions and thus clearly reduces 
flavor symmetry breaking. As such, we use hypercubic blocking to both estimate and 
reduce the effects of flavor symmetry breaking on our results. 

We generate our ensembles using the standard thin-link staggered fermion ac- 
tion (|5.47|) . We then apply hypercubic blocking to our ensembles and proceed with 
our analysis using both the thin-link and hypercubic-blocked ensembles. We have a 
greater trust in results arising from hypercubic-blocked ensembles, but by comparing 
those results to the thin-link results, we can estimate the magnitude of the systematic 
error that flavor symmetry breaking introduces. 

By using hypercubic blocking in this way, we are in effect using a different interac- 
tion matrix for our dynamical and valence quarks. Calculation of a partially quenched 
106 



quark-antiquark correlator after hypercubic blocking amounts to the expression: 

{,Qaai\nQbPj;m) 

= ^p'^qcd/ m] e-^^[^l detM(^,)[f/] M^^.^iV]-^,^^^ (5.117) 

While we admit that it is neither clear that this procedure has a clean continuum 
limit, nor that it corresponds to a well-defined field theory, it should at least prove 
useful in estimating the magnitude of fiavor-symmetry-breaking effects. 

It should be noted that, while we do not use hypercubic blocked links in our 
dynamical interaction matrix, an algorithm for doing so has since been proposed 



51, 52 . 



5.7 Sommer Scale 

When the creation of an ensemble begins, a set of values for the parameters in 
the lattice Lagrangian must be chosen. In our case these parameters are the gauge 



coupling, via P, and the unitless quark mass mg. As discussed in Section |5.1.3| , the 
lattice spacing a is not chosen directly, but rather is related to the other parameters by 
their renormalization group equations. In order to determine a, we use the ensemble 
to calculate some dimensionful quantity and then match the result to a physical 
measurement. 

A quantity which proves to be particularly useful for this purpose is the Som- 



mer scale |5J], also known as rg. Its popularity as a scale-setting quantity arises 
from the fact that it can be calculated on the lattice easily and with small statisti- 
cal error, its physical value is well determined, and its definition is equally valid in 
both quenched and unquenched simulations. No other quantity combines all of these 
advantages. 
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5.7.1 Static Quark Potential 

The definition of the Sommer scale is based on the static quark potential V{r), the 
potential energy function between two static quarks. The Sommer scale Tq is defined 
to be the separation distance r between two static quarks at which: 

r^F(r) = c = 1.65 (5.118) 

where F{r) is the force between the two quarks: 

F(r) = -^V{r) (5.119) 

The dimensionless parameter c = 1.65 is chosen so as to correspond to an inter-quark 
distance which is both convenient for calculation on the lattice and well explored 
phenomenologically. A second common choice is c = 1 which defines an alternative 
scale known as ri. 

The static quark potential dominates the internal physics of heavy quark mesons. 
Thus, its form can be gleaned from the spectra of both J/ip and T. All successful 
phenomenological potentials are in agreement as to the behavior of the static quark 
potential near the Sommer scale, giving: 

ro = 0.49fm (5.120) 

5.7.2 Wilson Loops 

On the lattice we calculate the static quark potential, and from it tq, using the 
expectation value of rectangular Wilson loops W^^*. In general a Wilson loop is any 
closed path of gauge-link matrices. In the case of rectangular Wilson loops, we define 
W^^.^ to be a loop with its lowest corner at site n and with a length of s gauge hnks 
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in the direction // and a breadth of t gauge hnks in the direction u: 



^v;n+sA+(^-l)!>^ ^M;n+(^-l)A^ (5.121) 

The action associated with a heavy point-hke color charge is simply the path- 
ordered integral of the gauge field along its trajectory. Prom this perspective the real 
trace of the Wilson loop, Re tr [1^.4^0*] , is thought of as the action associated with 
creating a quark-antiquark pair at the origin and instantaneously separating them a 
distance as, having the resulting states propagate an interval of time at, and finally 
instantaneously reuniting the pair. If we go to the limit of large time separations, 
the state which dominates the expectation value of the Wilson loop is two static 
quarks separated by a distance as. By its very definition, the static quark potential 
V{r — as) corresponds to the energy of this state: 

hm (Retr [W^*^*] > = ^e-"^^"^)* (5.122) 

where ^ is the amplitude for creating and annihilating the state. In practice, in order 
to maximize the information extracted from each configuration and thus minimize 
statistical error, we calculate the Wilson-loop expectation value by summing over all 
possible positions and allowed orientations of the loop: 

W^^'*^EE^M4;n (5.123) 
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5.7.3 Corrected Cornell Potential 

A standard ansatz for the continuum static quark potential is the Cornell poten- 
tial: 

aVc{as) — vo + vis + V2- (5.124) 

s 

This ansatz is based on two of the defining features of QCD: confinement and asymp- 
totic freedom. The second term is a string-tension term, which dominates at large s 
and causes the potential to be confining. The third term is a Coulomb term, which 
dominates at small s and results in asymptotic freedom. In our case we have ex- 
pressed the potential in terms of the unitless parameters Vi in order to facilitate our 
numerical analysis. 

The form for the continuum's Coulomb potential follows directly from the contin- 
uum's gluon propagator Gc{k): 

1 . r d'^k 



c 



- = 47ry -(^Gc{k) cossk, (5.125) 



where: 



Gc{k) = ^ (5.126) 
k ■ k 

and the momentum integration variable k is unitless. Recall that we have defined the 
displacement s to be in the direction /i. 

On the lattice, finite lattice spacing and the loss of rotational invariance alter the 
form of the gluon propagator at tree level, and thus modify the Coulomb potential. 
We denote the tree-level Coulomb potential for a system as: 

'■^ d^k 



-— Gx{k) COS sk, (5.127) 
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where X becomes L for the thin-hnk lattice and H for the hyper cubic-blocked lattice. 
The corresponding tree- level gluon propagators are |[55[]: 



GLik) 



p ■ p 



p ■ p 



where: 



Pi = 2 sin ■ 



hi 



(5.128) 



(5.129) 



The factor S^{k) accounts for the effects of the hypercubic blocking: 



where: 



(5.130) 



i^i = Pi 



1 + as (1 + as) - ^ (1 + 2a3) {p-p-p-) + ^ 11^3 



(5.131) 



and ctj are the hypercubic-blocking coefficients. 

The values used for the corrected Coulomb potentials are shown in Table p.l| . They 
were determined using a simple Monte Carlo integration of the momenta integrals 

(Em- 

We can now introduce a corrected potential: 



aVx(as) = vo + vis + V2- + V2 

s 



1 



SJ X s 



(5.132) 



where the parameter V2 allows for a correction of the continuum's Coulomb term. 

Once an ensemble's Wilson-loop expectation values have been calculated, we fit 
that two-dimensional data to the form: 



(Retr[Vr"^*]) = ^,exp< 



Vq + ViS + V2- + V2 

s 



1 

X s 



t 



(5.133) 



The range of s and t used in the fit must be chosen carefully in order that t remains 
large enough that the static-quark state dominates, while both s and t remain small 
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-1- 












.s. 


X 






C 


L 


H 


1 


1 


1.0817(12) 


0.7330(12) 


2 


0.5 


0.53793(73) 


0.48464(73) 


3 


0.33333 


0.34572(54) 


0.33359(53) 


4 


0.25 


0.25513(68) 


0.25187(68) 


5 


0.2 


0.20176(43) 


0.20071(43) 


6 


0.16667 


0.16868(71) 


0.16827(71) 


7 


0.14286 


0.14337(46) 


0.14319(46) 


8 


0.125 


0.12597(43) 


0.12586(43) 


9 


0.11111 


0.11157(38) 


0.11523(38) 


10 


0.1 


0.10009(46) 


0.10006(46) 


11 


0.09091 


0.09139(61) 


0.09137(61) 


12 


0.08333 


0.08355(66) 


0.08355(66) 



Tabic 5.1: Trcc-lcvcl Coulomb potential for the continuum C, thin-link lattice L, and 
hypercubic-blocked lattice H. The error listed is the statistical error of the Monte Carlo 
integration through which the values were determined. Because the same set of sample 
points was used for the Monte Carlo integration of the thin-link and hypercubic-blocked 
propagators, their error is highly correlated. 
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enough that the finite extent of the lattice does not come into play. The results of 
the fit are values for its free parameters: vq, vi, V2, V2, and an amplitude ^/g for each 
value of s included in the fit range. 

With these values in hand, the ensemble's lattice spacing can be determined via 
the Sommer scale. The Sommer scale tq is defined to be the distance at which the 
continuum static quark potential meets the criteria of ( |5.118|) . Thus after fitting, we 



drop the correction term from our potential determining Tq to be: 



Matching with phenomenological results ( |5.120|) , we find: 



'^o = aW (5.134) 



„-i = I 395 /i:5i±^ I MeV (5.136) 



or: 



Vl 

It is worth noting that, because the GL coefficients are unitless parameters, their 
calculation does not involve the lattice spacing directly. However, knowledge of the 
lattice spacing is still critical in that it allows access to such important values as the 
physical extent of an ensemble's lattice, the physical mass of the pseudo-Goldstone 
boson, and the expected strength of finite-lattice-spacing errors. 

5.8 Meson Propagators 

The lightest states produced by each quark bilinear are that bilinear's correspond- 
ing staggered meson and its parity partner. A two-point correlation function which 
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uses a bilinear as its creation and annihilation operator will thus, at large time separa- 
tions, be dominated by the propagator of these two states. Therefore, in the manner 
outlined in Section 14.31, the bilinear correlators allow us access to the masses of the 



staggered mesons, and in the case of the local pion, the decay constant as well. 

5.8.1 Bilinear Correlators 

In order to analyze meson propagators using lattice techniques, we must first 
express the bilinear correlators in a form which leaves them vulnerable to lattice 
investigation. 

We discuss here only local-bilinear correlators, 5 = JF, as they are significantly 
simpler than more general non-local correlators and are of primary importance in our 
study. For a similar discussion which encompasses non-local bilinears, see Appendix 



B.3. 



We begin with a simple correlator of the form presented in Section [4.3| . For 
our creation operator, or source, we use the bilinear with the quantum numbers of 
whatever meson we have chosen for study. For our annihilation operator, or sink, we 
use a wall of the same bilinear, summing over all positions on a time slice in order to 
restrict the annihilated states to those with zero spatial momentum: 



J2'^s,S:gJs,Sfl) (5.137) 



g 

gi=t 

The large time behavior of this correlator will be dominated by the propagator of the 
corresponding staggered meson and its parity partner. We denote the time separation 
of the correlator with the unitless parameter t, which takes on the integer values 
[0, U - 1] 
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We replace the single bilinear in our source with a linear combination of all local 
bilinears: 

^(E^^-^;^E^^.^;o) (5.138) 

Our source now overlaps with significantly more states, including all local mesons. 
However, with the bilinear form of the sink unchanged, only our desired states are 
annihilated, and thus the correlator continues to measure only the states we desire. 
Furthermore, we will find that this change allows us to calculate all local correlators 
simultaneously, minimizing our computational effort. 

For our source we switch from a single bilinear to a wall of bilinears at the ap- 
propriate time slice. This reduces the statistical error of our calculation without 
significantly increasing the computation time. The resulting factor of ^ = L1L2L3 is 



accounted for in Appendix |5.8.3 . 



Cs,s;t = ^(E^^.^;^ E (5.139) 

54=* h4=0 

The sum of all local bilinears has a simple form in terms of x ^^nd x, as demon- 
strated in Appendix |B.1| . Only a single contraction over the hypercube is non-zero, 
that at the lowest corner of the hypercube: 

J2^m = ^^XhXh (5.140) 
7^ 



Using ( |5.72|) and (|5.140| ), we put the correlator in terms of the fundamental lattice 



fermion degrees of freedom, x and x'- 



Cs,s;t = (E E(-)''^'"^"^^'+^^f+s E ^'^^h) (5-141) 



9 B Jl 
94=t /i4=0 
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The situation is complicated slightly by the fact that the fermion field at each 
site is not a scalar but rather a color vector. Within each bilinear these color vectors 
are contracted with an implied sum over color. That is, taking vn to be the fermion 

a 

degree of freedom with color a at site n: 



= (E E(-)''^'^^" E E^'^^v (5.142) 

a a ^ c c ' 

g ts a /i c 

94=* h4=0 



Via Wick contractions we compute the Grassmann integral over the fermion fields 
analytically, putting our correlator in terms of only gauge degrees of freedom. We do 
not allow the bilinears to self-contract by asserting that we only wish to investigate 
fiavor-non-singlet mesons. The contracted x ^i-nd x fields are transformed by the 
integral into the inverse of the fermion interaction matrix: 



Cs,S;t = (E E E(-)''^"'^'' E Xg+BXg+BXhXh) 

a,c g B ^ 



^ _ A) 

g4=t /i4=o 

E E E(-)''^'^^" E MmglBAM'[uf)~i^,, 

a,c g B J[ a,c 

a,c g B j; "'^ 

34=* /i4=0 

E E E(-)"^""""^^" E \M'[U]-Uh\') (5-143) 



a,c g B ^ 

94=* /l4=0 

1 '\ 1 In n o V~¥/~\^m T-\n-i- i m ^^-vry~ic^ ]\/fS\TT~\ ^ 



a,c 



where (M'^[f/] ) ^ has been put in terms of M'^[[/] as described in (|B.6|), and 5 is 



taken to be the binary four-vector with all elements set to one, such that 75 has the 
proper value under the definition ( ^.59| ). 

The bilinear correlator is now in terms of quantities which are calculable on the 
lattice: elements of the inverse of the fermion interaction matrix. Because we do not 
have access to the inverse interaction matrix itself, but rather use CG to apply the 
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inverse to a fermion field vector, calculating the individual elements of the inverse 
matrix, as appears to be required by the correlator, requires a prohibitive number 
NcV of applications of the inverse interaction matrix. 

In order to reduce the number of inversions required, we calculate the field vectors 

X{c}e. 

xL'^'^J2^'\.uCHvi (5.144) 



a,c 

h 



h4=0 

each of which requires only a single application of the inverse interaction matrix. 

T]^ is a set of noise field vectors for which each element equals a random unit-length 
phase. The introduction of this random-phase vector eliminates unwanted cross terms 
from the square of X^^^^. In the limit of an infinite number of noise vectors Nf 

J^^ Trl]^n^m = ^n,m (5.145) 

The sum for the off-diagonal elements is over an infinite number of random phases 
and washes out to zero. Thus: 

e i h "'^ / "'^ 

^14=0 /4=0 

= El^^[^]nir (5-146) 

n 

/l4=0 

The use of the random-phase vector has contracted M^[U]~'^ with itself, allowing us 
to avoid the calculation of its individual elements. 

In practice, we use only a single noise vector, — 1. This is allowed as the 
random phases between the cross terms within a single element of the sum over 
noise vectors is enough to wash out the position-off-diagonal terms in the square of 
X^^\ This single noise vector is effectively provided by the random phase which is 
automatically present at each lattice site due to the gauge hnk's local gauge freedom. 
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Our correlator becomes: 

^5,5;*= (EEE(-)"^^^'^^^^^El^i2Br) (5.147) 

c g B a " 

fl4=* 

where we have one field vector X^"^) for each color: 

h 

hi=0 

= (M^[t/]"V('=))n (5.148) 

a 

These Nc field vectors X'^'^'> each require an application of the inverse of the interaction 
matrix to calculate. To do so, we construct the field vector W^'^\ which equals one 
only at color c in the lowest corner of each hypercube on the time slice n4 = and 
zero elsewhere. The result of applying the inverse interaction matrix to W^'^^ is X^'^\ 

From inspection it is clear that the calculation of X*^'^-' is independent of both S 
and t. Thus, due to our use of the sum of all local bilinears as our source, we are able 
to calculate the correlator between all local bilinears at all time separations using 
only Nc inversions per configuration. 

In the case of the local pion, <S = 5 and '^'s+b,s+b;B = '■p'o,o;B — ^- The 

correlator becomes quite simple: 

c..« = (EEEEl-^SBr) (5.149) 

c g B a " 
34=* 

In practice we calculate the correlator multiple times, shifting our origin's time 
component and averaging over the results. For our study we repeated the calculation 
every four lattice steps in time. This allows us to extract the maximum amount of 
information from each configuration, minimizing our statistical error. 
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5.8.2 Meson Masses 



For large time separations the bilinear correlator will be dominated by the propa- 
gator of the lightest states it creates. In our case those are the appropriate staggered 
meson and its parity partner. For the moment we will ignore the parity partner. 

For a bilinear correlator with its source contained on the time slice = and its 
sink contained on the time slice = t, we see from ( [4. 181) that the large t behavior 
will be: 

i^e"''^* (5.150) 

where M is the mass of the staggered meson and is the amplitude for creating and 
annihilating the meson state. Note that we can define a unitless mass M = aM by 
absorbing into M the factor of a. Such unit-free quantities prove useful, as numerical 
procedures can only be used to determine dimensionless values. Throughout, we use 
the notation X to indicate some quantity X after the appropriate powers of lattice 
spacing have been absorbed so as to make X dimensionless. 

The lattice's periodic boundary conditions allow states to propagate from source 
to sink along both time directions. Accounting for this in our propagator, we find: 



As discussed in Section ^.8.1| , we do not use a single-time-slice sink, but rather 
a full-bilinear sink which lives both on time slice = t and = t + 1. Thus, in 
addition to this approximately doubling the magnitude of our correlator, the states 
must propagate different distances when being annihilated at the two time slices of 
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our sink: 



Accounting now for the negative-parity state, whose propagator oscillates in sign 
each time step, the final form for our large-time bilinear correlator is: 

Cs,T;t = ^+(1 + 6-^^^'+) (e-'^''^+* + e-'^^'^+(^^-i-*)) 

+ - e-^^-) (e-'^*^-* - e-«^^-{^4-i-*)) (5.153) 

where and are the amplitude for the creation and annihilation of the positive 
and negative parity states, and M+ and M_ are those states' masses. Recall that L4 
has been defined to be even. 

In the case of the local pion, there is no overlap with the negative parity state, 
^_ = 0: 

C5,5;t = <5 (1 + e-"*^-5) (e-'^^^-s* + e-'^^^-5(^4-i-t)) {^.IM) 

where is the amplitude for creating and annihilating the local pion state using 
our source and sink operators, and M^^ is the mass of the local pion. 

Thus, after completing a lattice calculation of 6*5^5;^, we can fit the results to the 
form described by ( |5.154| ) and extract the mass of the local pion. 



5.8.3 Meson Decay Constants 

While the local pion mass can be extracted without any consideration for the 
correlator's amplitude ^^s, the local pion decay constant /^^ can not. Thus, careful 
attention must be given to the normalization of the source and sink operators. 
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Recalling ( 4.18 ), we put our large-time correlator's amplitude in terms of our 



source and sink operator's overlap with the local pion state: 

C5,5;t = -i— (0|Osinkk5)(7r5|Osource|0)C:5,5;t (5.155) 

where [tts) is the dimensionless zero-momentum staggered local-pion state and: 

^5,5;t ={1 + e-*^-5) + e-*-5(^4-i-t)) ^^_-^^q^ 

From ( |5.141|) we identify our source and sink operators as: 

C^source ^ ^ XhXh ^ ^ XnXn 

jl it even 

/l4=0 714=0 

= ^nXnXn (5.157) 



n even 

714=0 



and: 



Osink = Y.Y1 i-T''-"'''Xg+BXg+B 

9 B 
9i=t B4=0 

= X] (5.158) 



n 
n4=t 



Our addition of an e„ factor to the source is allowed as the source is only non-zero 
at even sites. We consider only a single time slice, B4 = 0, of our two-time-slice sink 
operator, as the form of Cs.sjt already accounts for a second time slice. 

We can now put our correlator in terms of (0|e„x„x„|7r5), the overlap between the 
local pion state and a single-site operator with the local pion's quantum numbers. 
We find: 

V 

(OlCsourceks) = ^ (0 1 e„XnXn | TTs) (5.159) 

(OlOsinkks) = V(0|e„XnXn|vr5) (5.160) 
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and, thus: 



= T7^K0|enXnXnk5)pC:5,5;t (5.161) 

In the continuum, a meson's decay constant can be defined by the meson state's 
overlap with the appropriate weak interaction bihnear. In the case of the continuum 
pion, this definition is: 

^/2^M2 = {rriu + m^) {0\u-f5d\7i+) (5.162) 

where Ivr'*') is the continuum's zero-momentum tt'^ state. Recall that we are using the 
pion decay constant normalization in which /^^ ~ 92.4 MeV. In order to put our cor- 
relator in terms of the pion decay constant, we use the discretization correspondence 

(0|xZ75(i|7r+) — ^(0|e„XnXn|7r5) (5.163) 
a^v4 



Thus, combining ( |5.161|) , ( |5.162|) , and (|5.163| ), our correlator becomes: 



^5,5;t — Z 9 ^5,5;i 



fl Ml V 



^5,5;* (5.164) 



where f-,,^ is the local pion decay constant and f-,,^ is the corresponding dimensionless 
value. 

Knowing the normalization of our correlator, we now have the ability to extract 
both the local pion mass and decay constant. 

The local pion decay constant is protected by the staggered formulation's even- 
odd symmetry. This protection does not occur for any of the other fifteen staggered 
pions, as the staggered action breaks their corresponding axial vector symmetries. 
Thus, making a connection between the lattice propagator of one of these pions and 
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its continuum decay constant requires the calculation of a renormalization factor. 
The process is thus much more involved than for the local pion case. 
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CHAPTER 6 



PARTIALLY QUENCHED CHIRAL PERTURBATION 

THEORY 



In Chapter ^ we presented ChPT's prediction for the quark mass dependence 
of the chiral pseudo-Goldstone boson's mass and decay constant. In Chapter ^ we 
detailed lattice techniques for calculating that mass and decay constant. However, 
performing a set of lattice calculations using different values for the quark mass is 
extremely computationally cumbersome. As we are constrained by finite computer 
resources, we find it easier to vary only the valence quark mass. That is, we make 
use of the partially quenched approximation of Section ^.2.6 . 



ChPT is the low-energy effective field theory for the light bound states of un- 
quenched QCD. Thus, its predictions are not valid in the context of our partially 
quenched calculations. In order to bridge this gap between ChPT and partial quench- 
ing, we must make quantitative sense of the effects of partial quenching. The culmi- 
nation of this process is partially quenched Chiral Perturbation Theory (pqChPT), a 
low-energy effective field theory for the light bound states of pqQCD |58|, |59 |. 



Note that throughout this chapter, as in Chapter ^, we work in an infinite volume 
continuum. 
124 



6.1 Partially Quenched Quantum Chromodynamics 



As discussed in Section p.2.6| , pqQCD is not a well-behaved unitary field theory. 



Nonetheless, we are not prevented from proposing a quark content for pqQCD, nor 
from discussing pqQCD in a Lagrangian context. 

Partially quenched QCD clearly contains two types of fermionic quark flavors: 
dynamical quarks and valence quarks. As is illustrated by (|5.32| ), the dynamical 



quark field appears in the partition function's Boltzmann weight but does not appear 
in the operators whose expectation value we calculate. The valence quark field, on the 
other hand, appears only in these operators and does not contribute to the Boltzmann 
weight. The challenge before us is to construct a quark content for pqQCD in which 
these abnormal attributes natural consequence. 

Eliminating dynamical quarks from external operators is a simple matter of re- 
stricting ourselves to calculating only the expectation value of operators which involve 
valence quarks. By definition, this was our intent from the beginning. 

Eliminating the valence quarks' natural contribution to the theory's Boltzmann 
weight is more involved. We introduce a set of ghost quark flavors, scalar quark fields 
with incorrect spin statistics, with an interaction matrix equivalent to the valence 
quarks'. As the ghost quarks are bosons, their functional integral results in the inverse 
of the interaction- matrix determinant. If, for each valence quark flavor, we include 
a ghost quark flavor of equal mass, the two contributions to the theory's Boltzmann 
weight will cancel. The valence quark interaction-matrix determinant, which would 
naturally appear in the Boltzmann weight, is effectively eliminated. 

Thus, pqQCD can be thought of as having three types of quark flavors: dynam- 
ical quarks, valence quarks, and ghost quarks. The first two are standard fermionic 
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quarks and the last is bosonic. In our study, we use dynamical and valence quark 
flavors which are independently degenerate. The mass of the ghost quarks is then 
set equal to the valence quark mass. We represent the number of dynamical quarks 
as Nf and the number of valence quarks, and therefore ghost quarks as well, as A^^. 
From the perspective of this quark content, the pqQCD partition function, before the 
integration over quark degrees of freedom, is seen to be: 

ZpqQCD = J [VA,][Vq,][Dq,][Dqs][Vq,][Vq,][Vl] e- ^^"^--^^^ (6.1) 

where: 

■ZpqQCD = l^^'^[FfjLuF^''] 

+ q^ {YD^ + my) q^ + qs {^D^ + m,) + i]g {l^D^ + m^) qg (6.2) 

We use q^j to represent the vector of fermionic valence quark flavors with mass m^, 
qs to represent the vector of Nf fermionic dynamical quark flavors with mass m^, and 
(jg to represent the vector of bosonic ghost quark flavors with mass m^. Our use of 
rris here to denote the dynamical quark mass, and elsewhere to denote the Standard 
Model's strange quark mass, should remain clear via context. 

Performing the functional integral over quark degrees of freedom results in the 
cancellation discussed above: 

= J [DA^] e- det [^^D^ + m,] (6.3) 



We can see by comparison to (|5.32|) that the pqQCD partition function (|6.3|) has the 



appropriate form to correspond to the calculation of partially quenched expectation 

values. 
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The above construction of the quark content of pqQCD can also be understood in 
the context of perturbation theory. From this perspective, dynamical quarks appear 
only in closed quark loops, and not in the external legs of Feynman diagrams. Valence 
quarks, on the other hand, appear only in external legs and never in loops. To insure 
that dynamical quarks appear only in loops, we restrict ourselves to calculating matrix 
elements between valence-quark external states. To remove valence-quark loops, we 
again introduce ghost quarks. For every Feynman diagram which includes a valence- 
quark loop, there will be an identical diagram in which that loop is replaced by a 
ghost-quark loop. As bosons, the ghost-quark loops will appear with opposite sign. 
Thus, all diagrams which include valence-quark loops are canceled. 

It should be noted that, although we have conceived of a quark content and 
constructed a Lagrangian for pqQCD, the theory remains pathological. The advent 
of scalar quarks, with incorrect spin statistics, spoils unitarity. Furthermore, we 
are artificially restricting the physical states of the theory, removing dynamical- and 
ghost-quark external states. This effectively sets elements of the S matrix to zero, 
further ruining unitary. However, because it is the partially quenched approximation 
itself which is non-unitary, we can not hope to construct a unitary field theory with 
the necessary partition function. 

6.2 Partially Quenched Flavor Symmetry 

Just as the flavor symmetry group of QCD dictates the form of ChPT, so will the 
flavor symmetry group of pqQCD dictate the form of pqChPT. 

The flavor symmetry group of massless QCD, SU{Nf)L ® SU{Nf)R, follows di- 
rectly from its quark content. Massless pqQCD has an analogous flavor symmetry 
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group which also follows from the quark content proposed above. First, we collect 
the standard and ghost quarks into a single vector: 



1%. 



Q = [Qv Qs qg] 



(6.4) 



This mixture of fermions and bosons causes the flavor symmetry group to be a graded 
group, SU{N^+Nf\Nv)L^SU{Ny+Nf\Nv)R. In truth, there are subtleties which arise 
in the identification of pqQCD's flavor symmetry group. However, Sharpe demon- 



strates in that the above graded group generates the correct Ward identities for 
pqQCD. Thus, while it is not the exact flavor symmetry group of pqQCD, it is still 
appropriate to use in constructing pqChPT. 



6.3 Graded Groups 



Matrices which are elements of a graded group have the following distinguishing 
properties. The graded group element S G SU{N\M) can be written in block form: 



A C 
D B 



(6.5) 



where A and B are respectively N x N and M x M matrices of standard commuting 
numbers, while C and D are respectively N x M and M x N matrices of anticommut- 
ing numbers. In order to preserve the correct behavior of the adjoint, the complex 
conjugate of a product of these anticommuting numbers is defined to switch their 
order: 



(ab)* = b*a* 



(6.6) 
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where a and b are elements of C ot D. Finally, the super trace of a graded group 
element is defined as: 

sTrS = TrA-Tr5 (6.7) 
from which the definition of the super determinant is based: 

sdet S = exp<^ sTr In S \ = ^— ^ 6.8 

I -'J det B 

Using these properties for graded group elements, the graded group SU{Ny + 
Nf\Ny)L ® SU{Ny + Nf\Ny)fi has the appropriate behavior to act as the flavor sym- 
metry group for pqQCD's mixture of fermionic and bosonic quarks. 

6.4 Partially Quenched Chiral Lagrangian 

With the flavor sjnnmetry group for pqQCD in hand, we are now in a position 
to develop a low-energy effective field theory for the light bound states of pqQCD, 
mirroring closely the process in Chapter ^ 

We begin with an assumption that the full flavor symmetry of pqQCD is sponta- 
neously broken by a chiral condensate, from SU{Ny + NflNy)^ (g) SU{Ny + Nf\Ny)ji 
down to SU{Ny + Nf\Ny)v, just as occurs in full QCD as discussed in Section p. 1.4 . 



The result is set of Goldstone particles, pqQCD's light mesons. Note that we use the 
term mesons loosely, as several of these particles are fermions. The fermionic sub- 
set of Goldstone particles are those which correspond to broken symmetries whose 
generator's non-zero elements are contained by the sub-matrices C and D. 

Partially quenched ChPT will be a low-energy effective field theory for these Gold- 
stone particles. We collect the partially quenched mesons into a flavor-space matrix 
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multiplying each meson field tt" by its corresponding broken flavor symmetry gen- 
erator: 



(6.9) 



$ = TT V = 

_(kqq •rqq_ 

where (j)qq is an Ny + Nf x Ny + Nf matrix containing bosonic bound states of ordinary 
quark- antiquark pairs, is an N^ x N^ matrix containing bosonic bound states of 
ghost quark-antiquark pairs, and (p^^ and (pgg are respectively N^ + Nf x N^ and 
Ny X Ny + Nf matrices containing fermionic bound states of mixed quark-antiquark 
pairs. We further subdivide the matrix (f)gq, making its dynamical- and valence-quark 
blocks explicit: 



0, 



qq 



''vv Yvs 



(6.10) 



where the N^ x Ny matrix (p^^ holds the mesons which contain a valence quark- 
antiquark pair. As only valence-quark flelds are used in the operators whose expec- 
tation values we calculate, it is only this block of mesons whose masses we mea- 
sure directly in our partially quenched lattice calculations. Prom the perspective of 
pqChPT, all other mesons appear only at one-loop order and higher. From now on 
we will refer to a valence quark-antiquark meson as a pion, as it coincides with the 
state in our partially quenched staggered calculation to which we have designated 
that name. We deflne the unitary field matrix E as: 

E = e^'""""/^ e SU{N^ + Nf\N^) (6.11) 

The quark masses break pqQCD's flavor symmetry explicitly, granting the Gold- 
stone particles a light mass. We incorporate this symmetry breaking into pqChPT 
using a factor whose matrix structure is equivalent to the quark mass matrix: 

X = 2iiM = 2/xdiag(m^, . . .,m^, . . .,m„, . . .) (6-12) 

Nf 
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As stated above, we use A^^ valence quarks of mass ruy, Nj dynamical quarks of mass 
Tils, and Ny ghost quarks of mass m^. 

Considering all terms which respect both Lorentz and pqQCD's graded flavor 
symmetry, except via insertion of Ai, and expanding simultaneously in both me- 
son momentum and meson mass, we construct the NLO partially quenched chiral 
Euclidean Lagrangian: 

=5^Sc?PT = ^sTr[9,Sta'^S]-^sTr[Stx + xS] 
-Li(sTr[a^Sta'^E])' 

- L2 [sTt [S^S+a^S] sTr [S^S+S'^S] ) 

- L3 (sTr [a^s^a'^Ea^s^a^E] ) 

+ L4 (sTr [df^E^d^E] sTr [E+x + X^] ) 
+ L5(sTr[a^Et9''E(Etx + xS)]) 
-L6(sTr[Etx + xS])' 
-L7(sTr[Etx-xS])' 

-L8(sTr[EtxEtx] +sTr[xSxS]) (6.13) 

Note that this Lagrangian has the same form as that for standard ChPT ( p.24|) , with 
the only differences arising from our use of a graded flavor symmetry group. 

As long as we restrict our choice for the dynamical and valence quark masses 
such that they remain within the radius of convergence of pqChPT, the partially 
quenched chiral Lagrangian will accurately describe the two-dimensional quark mass 
dependence of the mesons' characteristics and behavior. 
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6.5 Next-to-Leading Order Meson Properties 

Using the NLO partially quenched chiral Lagrangian, we determine NLO expres- 
sions for the mass pT| : 

= zm^{4:iTff\ 1 + -^(2m„ - m^) In zm^ + -^{m^ - nis) 



+ zm^{2a^ — 0:5) + znigN f{2aQ — 0:4) j (6.14) 

and decay constant: 

fvv = /|l + —^{^v + rris) ln|(m^ + 771^) + zm^'^ + zm^A^/yj (6.15) 

of the partially quenched chiral pseudo-Goldstone boson containing two degenerate 
valence quarks, our pion. The quantity z is defined in (|]T0|). Note that in the 
unquenched case, = m.^ = rris, these expressions correspond exactly to those from 
ChPT: (|3l00D and (|3lOT| ). 



For the pqChPT expressions for the pion mass and decay constant in the case of 



arbitrary quark masses, the reader is referred to |]62[ . 

Note that in the limit of zero valence quark mass, rriy 0, the logarithmic term 
in M^^ diverges. Such divergent log terms in pqChPT expressions are often referred 
to as quenched chiral logs, or sometimes simply as quenched logs. They are directly 
analogous to the chiral logs of standard ChPT. The appearance of quenched chiral 
logs makes it clear that the true chiral limit of pqQCD is reached only by taking 
the dynamical and valence quark masses to zero simultaneously. In that case, the 
quenched log terms remain finite. 

Quenched logs also appear in expressions derived from the low-energy chiral theory 
of fully quenched QCD, quenched Chiral Perturbation Theory (qChPT). We will not 
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Figure 6.1: The quark-mass plane of pqQCD. 

discuss qChPT in detail, noting only that such logs are divergent when the only quark 
mass available for adjustment, the valence quark mass, is taken to zero. Thus, doubt 
is thrown onto the very existence of a chiral limit under the quenched approximation. 

6.6 Physical Results from Partially Quenched Calculations 

There are two critical points which allow pqChPT to act as a bridge between our 
unphysical partially quenched calculations and the physical GL coefficients. 

The ffist is that unquenched QCD is a subset of pqQCD. If we visualize all possible 
dynamical- and valence-quark-mass choices as defining the plane of pqQCD, as shown 
in Figure |6.1| , unquenched QCD is seen to lie on the diagonal line at which = m^. 
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The second point is that the two-dimensional quark mass dependence of pqChPT's 
Lagrangian is known and exphcitly included up to whatever order we have chosen to 
work. Thus, the GL coefficients are not functions of quark mass, but rather are 
constants throughout the quark-mass plane. Most importantly, they are constants as 
we move across the line of unquenched QCD. 

Thus, the GL coefficients of pqChPT and standard ChPT are the same. If we 
perform a partially quenched calculation and, via pqChPT, determine from the cal- 
culation values for the GL coefficients, we have in fact generated valid results for the 
GL coefficients of the physical world. 

Additionally, in a unquenched calculation, we are restricted to exploring quark 
mass dependence only along the unquenched line. In the case of the pion mass, it 
is clear from ( p.lOOl ) that such an exploration would grant us access only to the GL 



coefficient combination 2a8 — «5 + 2A^/a6 — A^/a4. On the other hand, the use of partial 
quenching grants us additional moment arms with which to extract information from 
our calculated quantities. In the case of the pion mass, varying the dynamical and 
valence quark masses independently allows for the determination of two coefficient 
combinations: 2a% — and 2aQ — a^. 

It should be noted that the dependence of the partially quenched chiral Lagrangian 
on the number of dynamical quark flavors Nf is not known. Thus, the GL coefficients 
are functions of Nf. As a consequence, in order for our partially quenched results to 
be considered valid results for the physical coefficients, we are forced to use a physical 
number of dynamical light quarks, Nf = 3, in our calculations. We in fact do just 
that. 
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In standard QCD, the can be thought of as obtaining its large mass via the sum- 
mation over quark-loop chains which appear in its propagator. From this perspective 
it is evident that the rj of fully quenched QCD will be lighter than expected, as any 
such quark loops have been removed. Thus, the rj must be included in any low-energy 
effective theory for qQCD such as qChPT. This inclusion introduces additional cou- 
plings to qChPT and destroys any correspondence of its couplings to the couplings 
of physical ChPT. While it was not initially known whether partially quenched QCD 



suffers from this same flaw, Sharpe demonstrates in |Q that the 77' of pqQCD is in 
fact heavy, and can be safely left out of pqChPT. Thus, the correspondence between 
the GL coefficients of pqChPT and ChPT is retained. 



6.7 Calculation of the Gasser-Leutwyler Coefficients 

A basic outline of our procedure for the calculation of the GL coefficients is as 
follows. We calculate the local pion's bilinear correlator, using N j = 3 staggered 
lattice techniques as detailed in Chapter |], for a range of dynamical and valence 
quark masses. We then fit those results, using the techniques described in Chapter |^ 
and over a fit range contained by the radius of convergence of pqChPT, to pqChPT's 
corresponding predictions for the correlator's quark-mass dependence. 

In order to facilitate a numerical fit of the lattice results to the predicted forms 
from pqChPT, those forms must be expressed in terms of dimensionless parameters. 
Making use of ( |5.164| ), we build pqChPT's prediction for the correlator's form: 



^5,5^, = ^^fY (1 + + e-*^-5(^^-i-*)) (6.16) 



5my 
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where, using ( |6.14| ): 



Ml^ = imy(47r/)^|l + —{2mv - ms) Inzmy + Jf^i'^v - ms) 

+ zmy {2as — as) + zmsNf (2aQ — a^j j> (6-17) 



using (|6.15| ): 



/vrs = /|l + —^{'^v + ms) \n^{mv + ms) + zmy'^ + imsA^/yj (6.18) 

and z and / are dimensionless parameters related to z and / by appropriate powers 
of the lattice spacing. The unitless dynamical and valence quark masses are denoted 
by ms and my. 

The products of this fit are values for the physical GL coefficients of ChPT. 
6.8 Constant Dynamical Quark Mass 

As discussed in Section [5.2.6|, changing the valence quark mass is significantly 



less computationally demanding than changing the dynamical quark mass. Thus, for 
those ensembles which have the smallest expected systematic error, and thus those 
which are the most computationally demanding, we use a single dynamical-quark- 
mass value, and vary only the valence quark mass. Fortunately, variation of the 
valence quark mass grants us access to the GL coefficient combination we are most 
interested in, the combination 2a^ — a^. 

The GL coefficient combination 2as — as originally became of interest because it is 



the unknown term in Am (|3.35| ), where Am is the difference in the NLO contributions 
to the physical pion and kaon masses ( p.34|) . Clearly, any such differences in the NLO 
terms arise due to the contrasting valence-quark content of the pion and kaon. Thus, 
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it is of no surprise that variation of only the valence quark mass grants access to the 
coefficient combination critical in the determination of Am- 

Without variation of the dynamical quark mass, the terms in ( |6.17| ) and ( |6.18| ) 



which depend on ms can not be accounted for. Thus, in cases where we have calcu- 
lated the local-bilinear correlator only along a line of constant dynamical quark mass, 
we absorb the unknown NLO ms dependence into the LO coefficients z and /. This 
change results in a deviation from the original expressions only when z and / appear 
in NLO terms, and thus the error due to this change occurs at NNLO. The resulting 
expression for the pion mass is: 

Ml^ = imy (47r/)2 1 1 + — (2my -m^) In zmy + Jf^i^v - ms) 

+ zniy {2as - a^) > (6.19) 



and for the pion decay constant is: 

/tts = •^j-^ + ^^{mv + ms) \n^{mv + ms) +imyy| (6.20) 
These forms are then used in (|6.16|) in place of ( |6.17|) and (|6.18| ). 
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CHAPTER 7 



DATA MODELING 

Throughout our study, in order to extract information from our lattice calcula- 
tions, we model the results to theoretical predictions. We present here details of that 
process. 

7.1 Lattice Measurements 

The product of a lattice calculation is the value of some operator evaluated under 
each gauge configuration of an ensemble. We refer to such an operator evaluation 
as a measurement, using the term in the loosest sense, as no actual experimental 
measurement is taking place. 

An operator's form will contain some number of adjustable parameters Xn, which 
we label as measure parameters. For each measurement a value is be chosen for each 
of these measure parameters. We codify the results of the measurements as: 

X2,X3, . . .) = {0{xi,X2,X3, . . .))[u]^ (7.1) 

where a denotes the configuration under which the operator was evaluated, and runs 
from 1 to N, where N is the number of configurations in the ensemble, and the 
various x„ designate the location in the operator's adjustable-parameter space at 
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which the measurement is made. Each measure parameter x„ takes on a discrete set 
of values where i enumerates the values for that measure parameter at which 
we have chosen to make measurements, running from 1 to NM;n- As an example, 
in the evaluation of a single bilinear correlator, there is one measure parameter, the 
bilinears' time separation, which takes on integer values in the range [0, L4 — 1]. 

We generally use the more compact notation Ua^x, where the index x enumerates 
all possible combinations of allowed values for the measure parameters. It runs from 
1 to Nm, where: 

NM = X{NM;n (7.2) 

n 

representing the total number of measurements made on each configuration. 
7.2 Data Blocking 



As a result of autocorrelation, as described in Section |4.2.4| , measurements on 
neighboring configurations of a Markov chain are correlated. However, to correctly 
account for the statistical error in our data, we require that each measurement be 
independent. In order to meet this requirement, we divide the Markov chain into 
blocks, each of which contains Nb = N/Nk neighboring configurations. We choose 
an Nb larger than the Markov chain's autocorrelation length. Thus, the average 
measurement within each block can be treated as a single uncorrelated measurement, 
independent of the results from other blocks. 

From now on we work only with this blocked data: 

yB;x = ^^ya;x (7.3) 

iVR 

It is worth noting that, for Gaussian distributed data, the standard error of our 
result is independent of the value chosen for A^^, assuming Nb is larger than the 
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chain's autocorrelation length but is still small enough that Nk remains large. In- 
creasing Nb reduces the variance between blocks but also reduces our sample size. 
For Gaussian distributed data, the two effects cancel exactly. 

7.3 Data Correlation 

We visualize all of the measurements on a single block as defining the location of 
a single point in a large- dimensional measurement space, where each possible combi- 
nation of measure-parameter values represents a dimension of the space. Prom this 
perspective, the index x of each measurement yB;x is seen as indexing the components 
of the point's position. 

The full set of measurements for an ensemble is seen as a cloud of points in this 
measurement space. We expect the data points to cluster around some central value 
and their density to die off as a Gaussian in all directions. The shape of this cloud, 
in addition to the rate of its falloff along the various directions, characterize the 
statistical error of our measurements. 

For each measurement there is an ensemble- average value: 



which expresses the width of the cloud along the measurement-space coordinate x. 
However, cr^ only effectively characterizes the measurements' error if the cloud is an 




(7.4) 



which corresponds to the expectation value of the operator O^- 



Additionally, each measurement has a variance: 




(7.5) 
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ellipse with axes that are flush with the coordinate directions of the measurement 
space. In such a case the width of the cloud along each coordinate direction is all 
that is needed to characterize the shape of the cloud. Yet, if the measurements are 
correlated, we will find that the axes of the eUipse do not fall along the coordinate 
directions. 

Invariably, measurements are correlated. That is, the measurements vary coher- 
ently between blocks. If, when a given measurement on a given block is found to 
be above the ensemble average, another nearby measurement tends to also be above 
the ensemble average, then those two measurements are said to be correlated. If the 
second measurement tends to be on the opposite side of the ensemble average as the 
first, the two measurements are anticorrelated. If the two measurements move above 
and below the ensemble average in an independent fashion, they are uncorrelated. 

The correlations between the Nm measurements are quantified by the covariance 
matrix C: 



Along the diagonal, the covariance matrix corresponds to the standard variance of 



between measurements. A large positive off-diagonal value corresponds to a strong 
correlation, while a large negative off-diagonal value corresponds to a strong anticor- 
relation. Completely uncorrelated data will have a diagonal covariance matrix. 

Often literature refers to a correlation matrix p rather than a covariance ma- 
trix. The correlation matrix is constructed to specify only the correlation between 




(7.6) 



a measurement, Cx,x = 



I, while the off-diagonal elements represent correlations 
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measurements, with all information concerning their error removed: 



1 



(7.7) 



1 1 two measurements are correlated, their value will tend to vary coherently. From 
the perspective of the data cloud, measurement points will tend to fall along a line 
diagonal to the coordinate directions of the two measurements. Thus, the data will 
form an ellipse with diagonal axes. 

The eigenvectors of the covariance matrix correspond to the off-diagonal axes of 
this ellipse. Additionally, the eigenvalues of the matrix correspond to the width of 
the cloud along those directions, effectively giving the data's standard error measured 
along the eUipse's axes. 

It is possible in certain pathological situations that data could form shapes more 
complex than a simple ellipse. The most common example is a boomerang-like shape, 
where the cloud may actually miss the ensemble average completely. Correct statis- 
tical error analysis of such data requires the use of moments higher than the second, 
and is beyond the scope of this discussion. 

7.4 Theoretical Models 

For each set of measurements taken, we have some theoretical form which we 
expect the data to match. Included in this theory function are some number of fit 
parameters q whose value we hope to glean from the data: 



/(ci,C2,C3, . . . ;a;i,X2,X3, . . .) 



(7.8) 



We also make use of the more compact notation jx{c\i C2, C3, . . .) or jxif)- 
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Given a set of values for the fit parameters, the theory function fx{c) specifies a 
point in measurement space. We would like to find the set of fit parameters which 
result in that point being as close as possible to the data's ensemble average. Addi- 
tionally, we would like to use a metric for this measure of closeness which accounts 
for the variation in magnitude of the data's statistical error with the direction of the 
distance in question. 

The correct metric for this purpose is the inverse of the covariance matrix. It 
weights distance in a given direction with a strength inversely proportional to the 
variance of the data along that direction. Using this metric the squared distance 
between our ensemble average and the theory function for a given set of fit parameters 
is: 



We are now left to minimize with respect to the fit parameters, and thus determine 
the optimal set of fit parameter values. 

By the nature of its definition, the covariance matrix is symmetric, and must be 
positive definite in the limit of a large number of configurations, Nk — ^ oo. Thus, 
Cholesky decomposition can be used to determine a triangular matrix L such 
that: 




(7.9) 



X 



y 



LL^ = C 



(7.10) 



or equivalently: 




(7.11) 
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By applying L ^ to our vector of correlated distances, we generate a vector of uncor- 
related distances weighted by the inverse of their uncertainty: 

r. = L-}y{yy-fy{c)) (7.12) 

In terms of Vx, is now simply: 

X^ = E^' (7-13) 

In cases where the number of data blocks Nk is not large relative to the number of 
measurements Nm, there is no guarantee that the covariance matrix will be positive 
definite. When it is not, we have no choice other than to set the matrix's off-diagonal 
elements to zero, effectively disregarding measurement correlation in the context of 
the fit. We find that this is often required in our study, as many of the fits have two, 
or even three, measurement parameters. Additionally, increasing our available data 
requires lengthening our Markov chains, a process which is exceedingly computation- 
ally intensive. These two factors often lead us to a very large Nm relative to our 
available Nk- 

Minimizing can be done using a nonlinear least-squares fitting technique. We 



begin each fit with the downhill simplex method pSl EM, also known as the amoeba 



method, and complete it using the Levenberg-Marquardt method [^, |65 
7.5 Jackknife Error 

We determine the statistical error of our fit parameters via a jackknife procedure. 
This involves repeating the entire fit process Nk times, each time leaving out a 
different block of data from the data set. The result is a vector of values q^-, where j 
runs from 1 to Nk, for each fit parameter q. The standard error for each fit parameter 
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is then: 



J 

where: 

When reporting the value of a fit parameter, we quote the statistical error using 
the form: 

(7.16) 

where q is the optimal fit parameter value found in the full fit and is given by 

When our desired result is some function of the fit parameters g[c\^ C2, C3, . . .), the 
variation in the fit parameters may be correlated such that the subsequent variation 
in the result is larger or smaller than what would be expected. In such a case we 
apply the jackknife procedure to the result itself: 



2 ^^K 



9 ^{9{cid,C2j,c^;j,...)- g) (7.17) 



J 

where: 

'9=^y2 ^2;j' C3;i, • • •) (7.18) 

3 

The result is then reported in the form: 

^(ci,C2,C3, . . .) ± (7.19) 
7.6 Specific Applications 

We applied the data modeling techniques described above in six aspects of our 
study. For clarity we present explicitly the measure parameters x„, fit parameters q, 
measurement operator Ox, and theory function /^(c) used in each fit. 
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Because we are only able to manipulate dimensionless quantities in our numerical 
fit procedure, measure and fit parameters must be unitless. Often we absorb into 
dimensionful quantities some number of powers of the lattice spacing a in order to 
construct a fit's required unitless quantities. We remind the reader of our use of X 
to denote the unitless version of some quantity X. 

7.6.1 Static Quark Potential 

We determine the static quark potential at a single spatial separation by fitting 
the rectangular Wilson-loop expectation value to an exponential: 

Xn-.t Ox : (El 

cr.s^.aV{as) Uc) : (pD 

7.6.2 Spatial Dependence of Static Quark Potential 

We extracted the form of the static quark potential by fitting the rectangular 
Wilson loop expectation value to its predicted form: 

ci : =c/s, fo, vx, V2, V2 fx{c) : (|5.133|) 

Note that a unique fit parameter ^ is used for each value of s in the fit range. 

7.6.3 Single Correlator 

In order to determine the mass and decay constant of the local pion for a single 
value of the valence quark mass, we fit the local-bilinear correlator to an exponential: 

Xn -.t a : (EHB) 

cr-L.M^, Lie) : i^m) 
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7.6.4 Quadratic Valence- Quark-Mass Dependence 

To determine the valence quark mass value at which the local pion mass equals 
the physical kaon mass, we fit the local-bilinear correlator at various valence quark 
masses to a phenomenologically motivated quadratic form for the local pion mass: 



Q : Mnv^ao, «i, ^2 /x(c) : (|9J), 

Note that an independent fit parameter sz/m^ is used for each value of my studied. 
For illustrative purposes, in one case we also fit the correlator to a cubic form. 



using the fit characteristics described above, but replacing ( |9.6|) with (|9.9| ). 

7.6.5 Chiral Valence- Quark-Mass Dependence 

To determine the GL coefficient combinations 2as — as and a^, we fit the local- 
bilinear correlator at various valence quark masses to the form for the local pion mass 
and decay constant predicted by pqChPT: 



Xn : t, mv : (|5A4^) 



Q : z, /, 2«8 - «5, «5 Uc) ■■ §M), O, (Pi 



For small volume ensembles, we are forced to add a constant term to the pion- 
mass form in the above fit. For these ensembles we use the fit characteristics described 
above, but replace (|6.19|) with ( p.lOp . 



7.6.6 Chiral Dynamical- and Valence- Quark-Mass Dependence 

In order to determine the GL coefficient combinations 2Q;g — a^, 2aQ — a^, a^, and 
04, we fit the local-bilinear correlator at various dynamical and valence quark masses 
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to the form for the local pion mass and decay constant predicted by pqChPT: 



In this case the measurement operator is not a function of one of the measure pa- 
rameters, the dynamical quark mass ms- Instead, ms is a parameter of the ensemble 
itself. Thus, measurements at different values of m5 come from different ensembles, 
and there can be no correlation between them. The elements of the covariance matrix 
which correspond to measurements at different dynamical quark masses are simply 
zero. 



Xn : t,mv,ms 
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CHAPTER 8 



ENSEMBLE DETAILS 



We present here the ensembles used in our study, reveahng the motivation behind 
and the parameters used in the generation of each ensemble. Table ^.l] displays a 
summary of this information. 

For all dynamical ensembles, a configuration was added to the ensemble every 
ten HMD trajectories along the Markov chain, with the first configuration in each 
ensemble being ten trajectories from the ensemble's starting condition. 

8.1 Primary Ensemble 

We refer to our primary ensemble as ensemble A, or after hypercubic blocking as 
ensemble A hyp. This ensemble has our largest lattice extent, 16^ x 32, and reasonable 
lattice spacing. As such, it is from A hyp that we will generate our quoted results. 
This is our only ensemble to include two Markov chains, one of which has an ordered 
starting condition, and the second of which has a disordered starting condition. Such 
a set of Markov chains is very helpful in determining the thermalization point, as in 
a sense the chains begin on opposite sides of the desired equilibrium ensemble. 
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Li 


L2 




Li 




ms 




start 


N 


Nt 


Nb 


Nk 


A 


16 


16 


16 


32 


5.3 


0.01 


3 





2250 


250 


200 


20 


















D 


2250 


250 






B 


12 


12 


16 


32 


5.3 


0.01 


3 





2250 


250 


200 


10 


C 


8 


8 


8 


32 


5.3 


0.01 


3 





10050 


250 


200 


49 


W 


8 


8 


8 


32 


5.115 


0.015 


3 





10300 


300 


100 


100 


X 


8 


8 


8 


32 


5.1235 


0.02 


3 





10300 


300 


100 


100 


Y 


8 


8 


8 


32 


5.132 


0.025 


3 


T 


10000 





100 


100 


Z 


8 


8 


8 


32 


5.151 


0.035 


3 


T 


10000 





100 


100 


Q 


16 


16 


16 


32 


5.8 







144 configurations 



Table 8.1: The lattice parameters for the ensembles used in our study. The starting 
conditions are denoted as O for an ordered start, D for a disordered start, and T for a 
thermal start. The values for Markov chain length N , thermalization point Nt-, and block 
length Nb arc given as trajectory counts. Nk corresponds to the number of blocks available 
in an ensemble. 



8.2 Finite- Volume Ensembles 



In order to study the magnitude of finite-volume error in our results, we use 
two ensembles with smaller lattice extent than our primary ensemble, leaving all 
other lattice parameters unchanged. We hope that results from the 12^ x 16 x 32 
ensemble, or ensemble B, deviate only slightly from those of our primary ensemble, 
demonstrating that the finite- volume error in ensemble A is under control. We expect 
that results from the 8^ x 32 ensemble, or ensemble C, will deviate significantly from 
our primary ensemble's results, as its physical volume is thought to be too small for 
our study. This should, however, make clear the effects of finite volume. Both of 
these finite- volume ensembles have ordered starting conditions. 
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8.3 Varying-Dynamical-Quark-Mass Ensembles 

To determine the GL coefficient combination 2aQ—a4, we require a set of ensembles 
between which the dynamical quark mass varies, but all other ensemble parameters 
are constant. In order to vary the dynamical quark mass, but leave the lattice spacing 
unchanged, we have made use of results from the Columbia group They have 



mapped out the Nf = 3, L4 = 4 finite temperature transition, determining several 



critical (3 and mg value pairs . We have generated four ensembles using these val- 
ues, and refer to them as ensembles W, X, Y, and Z. Ensembles W and X have ordered 
starting conditions, while ensembles Y and Z have thermalized starting conditions. 



using initial configurations from the ensembles of [37 . 

Because of the significant computational effort required to generate four distinct 
ensembles, we have used a small lattice extent. A reasonable physical volume is still 
obtained, as the lattice spacing corresponding to the ensemble parameters used is 
rather large. 

We also make use of ensemble W to study the magnitude of finite-lattice-spacing 
effects on our results. This ensemble has approximately twice the lattice spacing 
of our primary ensemble, yet has approximately equal physical volume. Thus, the 
deviation of its results from those of our primary ensemble give an indication of the 
effects of finite lattice spacing. 

8.4 Quenched Ensemble 

Quenched ChPT predicts a different functional form than pqChPT for the pion 



mass and decay constant's dependence on the valence quark mass [^]. Thus, it may 
be possible to distinguish a quenched ensemble from a partially quenched ensemble 
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via a study of that dependence. To this end, we include in our study a fully quenched 
ensemble with the same lattice extent as our primary ensemble and approximately 
equal lattice spacing. We refer to this ensemble as ensemble Q hyp, studying it only 
after hypercubic blocking. We do not make use of the predictions of qChPT, but 
instead treat the quenched ensemble as though it were partially quenched, with Nf — 
3 and mg = 0.01. If the effects of quenching are strong, the results from ensemble 
Q hyp should deviate significantly from those of ensemble A hyp. 
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CHAPTER 9 



ANALYSIS 



We present here the analysis of our ensembles and the processes which lead to our 
results. 

9.1 Thermalization and Block Length 

Figures p.l| through p.7| show the local pion's bilinear correlator C^^^-t calculated 



on the individual configurations of our ensembles. Such a plot makes obvious the 
effects of both thermalization and autocorrelation. For each ensemble we present the 
bilinear correlator at two time separations: t = and t = 15. 

The vertical dotted line in each plot demonstrates the value used for that en- 



semble's thermalization point Nt as described in Section [4.2.6| . The vertical range 



represents the block length Nb used for that ensemble as described in Section 17^2 . 

In the case of ensemble A, with its pair of Markov chains, we overlay the evolution 
of the two chains on a single plot. This makes the appropriate thermalization point 
quite clear. 
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9.2 Sommer Scale 



As described in Section p77[ the lattice spacing of each ensemble is set via the Som- 
mer scale. We present here our calculation of the lattice spacings, as well as a study 
of the effects of our fit-range choices on the results. We calculate the lattice spacing 
of the thin-link and hypercubic-blocked versions of each ensemble independently. 

In all value of Smin = 1 is used as the lower spatial bound of our fit ranges. 

9.2.1 Effective Potential 

In order to determine an appropriate minimum time separation tmin to use for our 
fit range, we calculate the effective potential Ves{s) for a range of time separations. 
The value of the effective potential between t and t + 1 is the exponential decay 
constant which describes the falloff of the Wilson-loop expectation value considering 
only those two time separations. Explicitly, the effective potential at t + | is defined 
as: 

(RetrTiy^^*]) ^ ^ 

aVos(s) 1 =ln7^ ^ ^ 9.1 

^ ' t+\ (Retr[W^^x*+i]) ^ ' 

The effective potential is expected to vary rapidly for small t, where the asymptotic 
static quark state is infected by higher-energy states. Beyond some larger t the value 
is expected to plateau at the static quark potential. We chose tmin such that our fit 
range spans only the region in which the static quark potential is uncontaminated. 



Figures p.8| through |D.15| show the effective potential over a range of time sep- 
arations t and at two spatial separations s for each ensemble. The chosen spatial 
separations bracket, or nearly bracket, that ensemble's result for r^/a. The dotted 
vertical line in each plot represents the value of tmin chosen for the final fit. The error 
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bars displayed are the result of a jackknife analysis. Those points whose error bars 
span their plot's vertical range have been dropped. 

For ensembles W through Z, and their hypercubic blocked counterparts, our choice 
of tmin = 2 may not seem optimal in light of the effective potential plots. However, 
tmin must bc low cuough so that the statistical errors of the static quark potential are 
under control out to a spatial separation of s = 4. This requirement arises from the 
fact that we are using a four-parameter form for the s dependence of the Wilson-loop 
expectation value ( |5.132| ). Thus, without four well-determined values for the static 
quark potential along the s direction, the fit is under-determined. Choosing t^i^ = 3, 
which would perhaps seem advisable based on the effective potential's behavior, causes 
the statistical error of the static quark potential at s = 4 to be very large. Thus, we 



choose tmin = 2. In Section |9.2.3| we present the dependence of the resulting value of 
To/a on our choice of tmin- This dependence demonstrates that our choice of tmin = 2 
for ensembles W through Z in not unreasonable. 

If we were not using a corrected ansatz for the static quark potential, an under- 
determined spatial dependence would not be as significant an issue. When calculating 
jackknife error for such an under-determined potential, the results for the fit param- 
eters will vary wildly. However, as long as the statistical error in the vicinity of the 
Sommer scale is under control, these variations will be correlated such that the cor- 
responding variation in rg/a will remain small. For our corrected potential, we do 
not have this luxury. In this case, we ignore a term in the static quark potential's fit 
form when determining the Sommer scale. It is critical that the coefficient for that 
term V2 be well determined by the fit, as the correlation of its error with that of the 
other fit parameters is discarded, and thus is not available to stabilize ro/a. 
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It follows from the definition of the rectangular Wilson loop that its expectation 
value has a symmetric dependence on s and t. Based on our ansatz for the static 
quark potential ( p.l32| ), we can only expect the expectation value to fall off in time 
as a single exponential for time separations t larger than some value Sstr, where Sstr 
is the spatial separation beyond which the static quark potential is dominated by its 
string-tension term. The Sommer scale tq is specifically defined so that it falls in the 
transition region between a Coulomb-like and a string-like quark potential. Thus, it 
is clear that one should choose a value for tmin which is just greater than ro/a. In 
all cases we choose tmin to be just greater than the value determined for rg/a on the 
corresponding hypercubic blocked ensemble. 

9.2.2 Static Quark Potential 



For each ensemble the Wilson-loop expectation value is fit to our ansatz ( p.l33| ) 



over the range of spatial and time separations bound by Smin, Smax, ^min, and 
The value of ro/a is then determined via ( p. 134 ). 



In Figures D.16 through D.30 we display the results of these fits, plotting the 



determined corrected static quark potential, defined as: 



aVcorrias) = Vq + ViS + V2- (9.2) 

s 



The parameters Vq, Vi, and V2 are determined by the fit. The parameter V2 is also 
determined by the fit, but its term is dropped when calculating the corrected static 
quark potential. 

Included in the plots are the results of independent fits of the Wilson-loop expec- 
tation value along lines of constant s. The Wilson- loop expectation value is fit to the 
156 



form: 



(9.3) 



for a single spatial separation s and over the range of time separations t from imin to 
^max- The result is a value for the static quark potential at a single spatial separation 
aVias). These values are represented on the plots by x's. The plots' diamonds 
correspond to the corrected static quark potential, which in this case is defined as: 



where V2 is determined by the full fit as discussed above. 

It should be stressed that the corrected static-quark-potential curve is not the 
result of a fit to the corrected static-quark-potential points. The curve and points 
are only related in that they are both the result of fits to the same expectation- 
value data. Agreement between the corrected static-quark-potential curve and points 
demonstrates that our ansatz for the s dependence of the Wilson-loop expectation 
value is appropriate. 

The error bars are determined via a jackknife analysis. For clarity, error bars are 
only shown on the corrected potential points, although equivalent error bars on the 
corresponding uncorrected potential points would be appropriate. Each fit curve's 
one-sigma range is bound by dotted fines. This range is the result of a jackknife 
analysis of the curve's value at each point along the horizontal axis. Due to the 
large number of Wilson-loop expectation values used for the fit, the data's covariance 
matrix fails to be positive definite for all ensembles. Thus, in all cases we use a 
diagonal covariance matrix. 




(9.4) 
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On each plot the vertical dashed line corresponds to the ensemble's determined 
value of ro/a, while the shaded region corresponds that result's jackknife error bars. 



9.2.3 Dependence on tmin 

In order to directly study the tmin dependence of our results, we repeat the fits 
using a range of values for tmin- AH other fit range parameters are left unchanged. 
The fruits of this analysis are presented in Figures p.31| through |D.34|. In each plot 



the train uscd in the final fit is denoted by a filled diamond, which thus corresponds 
to that ensemble's determined value for ro/a. For several ensembles, especially the 
thin-link ensembles, very few points are given. For those values of tmin at which a 
point is not given, the fit either fails entirely, results in an imaginary value for ro/a, 
or produces statistical error bars which span the plot's vertical axis. 



9.2.4 Dependence on t 



max 



In order to directly study the tmax dependence of our results, we repeat the fits 
using a range of values for tmax- AH other fit range parameters are left unchanged. 
The products of this analysis are presented in Figures p.35| through p.38|. In each plot 



the traax uscd in the final fit is denoted by a filled diamond, which thus corresponds to 
that ensemble's determined value for ro/a. Due to increasingly large statistical error 
in the Wilson-loop expectation value for large time separations, the fit results have 
a weak dependence on tmax- Ensemble C hyp shows the strongest dependence, most 
likely due to significant finite-volume effects. 
158 







^min 


^max 


^^0 


Vl 


V2 


V2 


A hyp 


7 


4 


10 


0.243(16) 


0.1031(31) 


-0.336(21) 


-0.352(31) 


B hyp 


5 


4 


10 


0.293(24) 


0.0924(41) 


-0.396(33) 


-0.432(48) 


C hyp 


4 


4 


10 


0.374(15) 


0.0649(33) 


-0.479(19) 


-0.523(26) 


A 


7 


4 


10 


0.77(13) 


0.097(45) 


-0.16(11) 


-2.4(24) 


B 


5 


4 


10 


0.81(16) 


0.081(55) 


-0.30(17) 


-1.0(33) 


C 


4 


4 


10 


0.81(10) 


0.074(36) 


-0.33(13) 


-0.6(23) 


W hyp 


4 


2 


10 


-0.239(50) 


0.4715(99) 


-0.026(64) 


-0.119(89) 


X hyp 


4 


2 


10 


-0.129(40) 


0.4393(76) 


-0.152(52) 


-0.278(74) 


Yhyp 


4 


2 


10 


-0.144(47) 


0.4322(90) 


-0.114(61) 


-0.205(86) 


Z hyp 


4 


2 


10 


-0.128(43) 


0.4190(79) 


-0.129(57) 


-0.233(82) 


W 


4 


2 


10 


0.65(21) 


0.384(72) 


0.05(18) 


-4.9(38) 


X 


4 


2 


10 


0.33(19) 


0.487(67) 


-0.10(17) 


-0.4(36) 


Y 


4 


2 


10 


-0.05(20) 


0.613(69) 


-0.48(17) 


7.4(36) 


Z 


4 


2 


10 


0.11(20) 


0.548(72) 


-0.39(16) 


5.1(37) 


Q hyp 


7 


4 


10 


0.169(12) 


0.1043(22) 


-0.261(15) 


-0.290(22) 



Table 9.1: Fit-range limits and the results for a subset of the fit parameters from the 
corrected static-quark-potential fits. 



9.2.5 Dependence on s 



max 



In order to directly study the Smax dependence of our results, we repeat the fits 
using a range of values for Smax- All other fit range parameters are left unchanged. 
The outcome of this analysis is presented in Figures p.39| through p.42| . In each plot 
the Smax used in the final fit is denoted by a filled diamond, which thus corresponds to 
that ensemble's determined value for ro/a. Due to increasingly large statistical error 
in the Wilson-loop expectation value for large spatial separations, the fit results have 
a weak dependence on Smax- Again, ensemble C hyp shows the strongest dependence. 

159 



9.2.6 Results 



Tables |9.1| and |9.2| compile the results of the corrected static-quark-potential fits. 



Table presents the fit ranges used, and displays each ensemble's determined values 
for a relevant subset of the fit parameters. Table |9]^ presents the corresponding values 
for r^/a and the inverse lattice spacing. The quoted uncertainties are from a jackknife 
analysis of the statistical error. 

In order to gauge the effects of using a corrected potential, we repeat the fits using 
an uncorrected static quark potential, fixing the correction term's parameter to zero. 



^2 = 0. The results of these uncorrected fits are given in Table 9.2 



Note that between the thin-link ensembles W through Z, the results of the cor- 
rected fit show uncontrolled variation. This is most obvious in the values obtained for 



the correction term's parameter V2, and can be seen clearly in Figures p.26| through 
D.29| and in Table |9.1| . This is the result of a poorly determined s dependence for the 
static quark potential, resulting in an extremely fiat minimum for along a path 
through fit-parameter space. This leads to large, but highly correlated, variations in 
the fit parameters. Because we ignore the correction term when determining rg/a, 
these correlated errors do not have a chance to counteract one another. The end re- 
sult for ro/a is large statistical errors within each ensemble and large variations in the 
value between ensembles. For the corresponding uncorrected static-quark-potential 
fits, the situation is improved twofold. The ansatz for the potential contains fewer 
terms, and thus the resulting fit parameters are better determined. Additionally, 
all terms are retained when determining ro/a, and thus any correlated error has an 
opportunity to cancel in the final result for ro/a. This improved situation is clear 
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ro/a 

uncorr corr 


a^^ (MeV) 
uncorr corr 


A hyp 
B hyp 
C hyp 


3.507(23) 3.570(27) 
3.581(68) 3.684(44) 
3.757(28) 4.248(78) 


1384.1(90) 1409(10) 
1414(27) 1454(17) 
1483(11) 1677(31) 


A 
B 
C 


3.17(12) 3.9(10) 
3.69(24) 4.1(17) 
3.94(29) 4.2(12) 


1251(49) 1550(400) 
1457(95) 1610(670) 
1560(110) 1670(480) 


W hyp 
X hyp 
Y hyp 
Z hyp 


1.8809(37) 1.856(19) 
1.9076(33) 1.847(17) 
1.9293(35) 1.885(19) 
1.9547(40) 1.905(19) 


742.4(15) 732.5(73) 
752.9(13) 728.9(69) 
761.5(14) 744.2(75) 
771.5(16) 751.9(75) 


W 
X 
Y 
Z 


1.766(12) 2.10(30) 
1.763(11) 1.78(21) 
1.7809(93) 1.38(17) 
1.812(10) 1.51(19) 


697.2(47) 830(120) 
696.0(44) 704(84) 
702.9(37) 545(68) 
715.1(41) 598(76) 


Q hyp 


3.570(12) 3.650(19) 


1409.2(47) 1440.5(76) 



Table 9.2: The final results for the Sommer scale and lattice spacing from the corrected 
and uncorrected static-quark-potential fits. 



from the relative consistency between the uncorrected results for the Sommer scale 
of ensembles W through Z, as shown in Table |9.2| . 

One advantage of hypercubic blocking, a reduction in statistical error, is evident in 
our results for the Sommer scale. As such, we have a greater trust in our hypercubic- 
blocked results. In following sections we use an ensemble's hypercubic-blocked lattice 
spacing for both the thin-link and hypercubic-blocked versions of the ensemble. 

We suspect the variation in the determined lattice spacing for ensembles A through 
C is due not to true changes in the lattice spacing, but rather to our limited ability 
to determine the lattice spacing on the smaller-volume ensembles. Thus, in order to 
clarify the effects of finite volume in other quantities, we use, in the following sections, 
the lattice spacing determined for ensemble A hyp for ensembles B and C as well. 
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9.3 Effective Mass 



We can extract the value of the local pion mass and decay constant from the 
bilinear correlator C^^^-t because, for large time separations t, it contains only a single 
state, the static local pion. Thus, we must determine some time separation tram 
beyond which we will assume the pion state is uncontaminated. 

In order to make this choice for each ensemble, we calculate the effective mass 
of the local pion over a range of time separations. The effective mass Meg at t + | 
is the mass which describes the exponential falloff of the correlator, considering only 
the correlator's value at the time separations t and t + 1. In practice we determine 
the effective mass by fitting the correlator at a single valence quark mass to the form 
in (|5.164|) , using the value of the correlator only at four time separations: t, t + 1, 



L4 — 2 — t, and L4 — 1 — t. 



The resulting values for the effective mass are presented in Figures p.43| through 



D.46| . The error bars shown are the result of a jackknife analysis. In each case the 
effective mass is high for small t, where the correlator still contains higher-energy 
states, and then plateaus for large t. The value of tmin chosen for each ensemble is 
represented in the plots by a vertical dotted line. Those values can also be found in 



Table 9.3. All other fits of the correlator in this study use these values for tmm, fitting 



the correlator only over the range from tmin to L4 — 1 — tmin- 

9.4 Pion Mass 

Once we have determined an appropriate tmin, the local pion mass at a given 
valence quark mass can be ascertained by fitting the bilinear correlator C^^^-t to its 
expected exponential falloff (|5.164| ) in t. Table |9.3| shows the results of these fits 
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my = 


0.01 


my = 


= ms 




^min 




(MeV) 




M^5 (MeV) 


A hyp 


7 


0.19769(57) 


278.5(21) 






B hyp 


7 


0.20235(93) 


285.1(24) 






C hyp 


7 


0.335(11) 


471(16) 






A 


7 


0.29704(31) 


418.5(30) 






B 


7 


0.29792(37) 


419.8(30) 






C 


7 


0.3483(32) 


490.8(57) 






W hyp 


4 


0.25189(19) 


184.5(18) 


0.30659(19) 


224.6(22) 


X hyp 


4 


0.25143(25) 


183.3(17) 


0.35134(24) 


256.1(24) 


Yhyp 


4 


0.25098(23) 


186.8(19) 


0.39029(25) 


290.5(29) 


Z hyp 


4 


0.25091(26) 


188.7(19) 


0.45804(26) 


344.4(34) 


W 


4 


0.258565(96) 


189.4(19) 


0.315645(98) 


231.2(23) 


X 


4 


0.25955(11) 


189.2(18) 


0.36470(11) 


265.8(25) 


Y 


4 


0.260230(99) 


193.7(20) 


0.40742(11) 


303.2(31) 


Z 


4 


0.26140(12) 


196.5(20) 


0.48119(12) 


361.8(36) 


Q hyp 


7 


0.19472(70) 


280.5(18) 





Table 9.3: Minimum time separation and local pion mass at my = 0.01 and my = ms- 
For those ensembles in which ms = 0.01, the results are not repeated. 



using the valence-quark-mass values my = 0.01 and my = ms- For those ensembles 
in which ms = 0.01, the results are not repeated. The uncertainties quoted are the 
result of a jackknife analysis of the statistical error. In the case of the dimensionful 
result, the statistical error of the inverse lattice spacing is added in quadrature. We 
use a diagonal correlation matrix for these fits. 

Ensembles A and B generate very similar pion masses, indicating that the local 
pion mass is relatively free of finite- volume effects. The pion mass of ensemble C is sig- 
nificantly larger. This indicates that the small volume of the ensemble is constricting 
the pion, preventing it from relaxing to its lowest energy state. 

As discussed in Section p.4| , pqChPT is only valid for quark masses within some 
radius of convergence from the chiral limit. In order to gain some insight into the 
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physical magnitude of our quark masses, we compare our calculated pion masses to 
the experimental mass of the the Standard Model's light mesons. The experimental 
mass of the Standard Model's pion is M^o = 135.0 MeV, while the kaon mass is 
Mk+ = 493.7 MeV [^. From Table P73| , it is clear that the quark masses we use are 
larger than those of the up and down quarks. Yet, the ensembles' pion masses remain 
below the kaon mass. In the case of ensemble A hyp, the pion mass is well below 
the kaon mass. This leaves us hopeful that pqChPT is valid within the quark-mass 
range of our study. It should be noted that our coarse lattice spacings, especially in 
ensembles W through Z, are a large part of why our pion masses remain low. 
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9.5 Kaon Quark-Mass Threshold 

For staggered lattice fermions the quark mass is multiplicatively renormalized. 
While we do not wish to calculate the renormalization factors exactly, we do require 
some method for comparing the quark mass between ensembles. Thus, for each 
ensemble we calculate what we identify as the kaon quark- mass threshold rriQ^, the 
valence quark mass at which the local pion mass equals the Standard Model's kaon 
mass. At that point, the physical valence quark mass should be on the order of half 
the strange quark mass. 

In order to determine mg^, we require an ansatz for the pion mass's valence- 
quark-mass dependence which is accurate up to large valence quark mass. While 
pqChPT predicts a form for the pion mass, it is not appropriate to use it in this case, 
as the fit includes valence quark masses which are beyond the expected radius of 
convergence of pqChPT. Instead, we simply fit the squared pion mass to a quadratic 
polynomial, a form which it follows very well up to the largest valence quark mass we 
investigate. 

The phenomenological fit form we use for this purpose alone is: 

C,,,.,t = ^my (e-^-«* + e-*-5(^4-i-t)) (9.5) 

where: 

= oo + aiuiv + a2my (9.6) 

and s^my is a set of independent fit parameters, one for each valence-quark-mass value 
studied. We should stress that this is clearly not the form predicted by pqChPT. We 
use it here only because we require a good fit of the local pion mass up to large valence 
quark mass. The results of these fits are used for nothing other than an interpolation 
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CbQ dl (I2 


rriQK 


A hyp 
B hyp 
C hyp 


0.00223(11) 3.684(20) -1.99(22) 
0.00401(29) 3.694(37) -2.15(54) 
0.0671(76) 4.59(21) -10.1(25) 


0.03331(13) 
0.03276(20) 
0.0125(17) 


A 
B 
C 


0.00435(18) 8.558(15) -21.99(24) 
0.00503(20) 8.547(15) -21.85(24) 
0.0359(25) 8.744(49) -27.49(52) 


0.014365(21) 
0.014295(27) 
0.01027(27) 


W hyp 
X hyp 
Y hyp 
Z hyp 


0.00378(12) 6.0518(72) -2.778(47) 
0.00395(13) 6.0084(71) -2.654(45) 
0.00406(12) 5.9758(73) -2.617(49) 
0.00434(15) 5.9523(81) -2.669(51) 


0.077151(55) 
0.078381(63) 
0.075492(63) 
0.074160(64) 


W 
X 
Y 
Z 


0.001263(59) 6.6238(40) -4.911(28) 
0.001388(55) 6.6679(40) -5.100(30) 
0.001509(56) 6.6935(40) -5.204(28) 
0.001612(72) 6.7502(43) -5.456(29) 


0.072242(28) 
0.072598(25) 
0.069284(26) 
0.067282(27) 


Q hyp 


0.00269(22) 3.524(21) -2.25(27) 


0.03327(12) 



Table 9.4: Kaon quark-mass threshold and results for a subset of the fit parameters from 
the quadratic fits of the pion mass's valence-quark-mass dependence. 



to determine the point at which the local pion mass crosses the Standard Model's 
physical kaon mass: 



-ai + Jaj - 4a2(ao - (ax 493.7 Me V)^) 
= ^ (9.T) 

where a is the lattice spacing used for the ensemble. 

Figures p. 47| through p.50| display, and Table |9.4] summarizes, the results of these 
fits. The curve in each plot demonstrates the valence-quark-mass dependence deter- 
mined by the fit, and is generated by using (|9.6| ) and the resultant values for the 
fit parameters. The diamonds are the result of independent fits of the correlator at 



separate valence quark masses, using the same method as was used in Section |9.4| . It 
should be stressed that the full quadratic fit of the valence-quark-mass dependence 
is not a fit to these points. The points are related to the fit only in that they are 
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both derived from the same correlator data. Agreement between the fit curve and the 
diamonds demonstrates that our phenomenological ansatz for the local pion mass's 
valence-quark-mass dependence is appropriate. The valence-quark-mass values at 
which pion-mass points are given correspond to the values used in the quadratic fit. 

The error bars were determined via a jackknife analysis. Each fit curve's one-sigma 
range is bound by dotted lines. This range is the result of a jackknife analysis of the 
curve's value at each point along the horizontal axis. In most cases this region is 
thinner than the curve's line, and thus is not easily visible. Uncertainty in the lattice 
spacing is not taken into account. Instead, the plots' vertical axes are simply rescaled 
using the lattice spacing's central value. We use a diagonal correlation matrix in the 
fit, as the full correlation matrix is not positive definite for several of the ensembles. 

The dashed vertical line appearing in the plots corresponds to that ensemble's 
determined value for mg^ , while the shaded region corresponds to the uncertainty in 
that result, as determined by a jackknife analysis of the statistical error. In many cases 
this region is thiner than the line itself, and is not readily visible. This error, which is 



repeated in Table does not take into account uncertainty in the ensembles' lattice 
spacings. 

Because finite-volume effects inflate the pion mass in ensembles C hyp and C, the 
kaon quark-mass threshold for these ensembles is unexpectedly small. Thus, in order 
to clarify finite- volume effects in other quantities, we use the value of determined 
for A hyp for both ensembles B hyp and C hyp. Similarly, A's value for mq^^ is used 
for ensembles B and C. 
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9.6 2^8 — 



In order to determine the value of the GL coefficient combination 2as — 05, we fit 
the correlator C^^^-t over a range of time separations t and valence quark masses my 
to the form predicted by pqChPT: ( |6.16| ), (|6.19| ), and ( |6.2(J| ). The result is values for 



the fit's free parameters, one of which is 2as — 05. 

9.6.1 Valence- Quark-Mass Cutoff 

pqChPT is expected to accurately model the behavior of the local pion only within 
some radius of convergence of the chiral limit. Thus, we should only fit the correlator 
to pqChPT's predictions at and below some cutoff in the valence quark mass Amy 
While we expect the appropriate cutoff to be below the kaon quark-mass threshold, 
there is no clear a priori value. As such, we turn to the data to determine our cutoffi 

Figure |D.51| plots degree of freedom for fits of ensemble A hyp's correlator 



data to the predictions of pqChPT against various choices for the cutoff Amy A 
diagonal correlation matrix was used in these fits. The error bars are the result of a 
jackknife analysis. We choose as our value for Amy the cutoff at which per degree 
of freedom is closest to one: A^y = 0.025. In Figure p.51|the chosen cutoff is denoted 



by a filled diamond. This cutoff is then used for ensembles A hyp through C hyp, as 
well as for ensemble Q hyp. 

We know from our study of the kaon quark-mass threshold mq^ that the quark- 
mass renormalization factor of the thin-link ensembles is significantly different than 
that of the hypercubic-blocked ensembles. Thus, it would not be appropriate to use 
the same valence-quark-mass cutoff for the thin-link ensembles. Instead we determine 
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a cutoff which is consistent between the thin-hnk and hypercubic-blocked ensembles 
using tuq^. 

The ratio of the valence- quark- mass cutoff over the kaon quark-mass threshold 
for ensemble A hyp is hmv/^QK ~ 0-751. We choose the cutoff for ensemble A to be 
that which gives the closest ratio beyond this value. We choose Kmv = 0.0125, which 
corresponds to a ratio of K^v/rnQj^ = 0.870. This cutoff is then used for ensembles 
A through C. We select the cutoff for ensembles W hyp through Z hyp similarly, 
choosing Amy = 0.06. Between the four ensembles, this gives an average ratio of 
^mv/^QK — 0.0763. For ensembles W through Z, we use the cutoff Amy = 0.055, 
which between the four ensembles gives an average ratio of Amy/mQj^ = 0.782. 

Given the evident correlation present in the data, and that our calculations of 
do not take that correlation into account, it could be argued that, while clearly 
minimization of this is appropriate, its actual value is meaningless. As such, there 
is no compelling reason to believe that a choice of valence-quark-mass cutoff which 
results in a P^r degree of freedom of one is appropriate. However, as we are left 
with no other quantitative evidence as to an appropriate cutoff choice, we feel that 
it is important to at least use a consistent and systematic method for resolving the 
choice. 



In section |9.6.7| we study the effect of these cutoff choices on the resulting value 



of 2as — as and find that it is a significant source of systematic uncertainty. 



9.6.2 Pion Mass and Decay Constant 



Figures p.52| and |D.53| present results from the fit of ensemble A hyp's correlator to 



the predictions of pqChPT. These plots of the dependence of the pion mass and decay 
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constant on the valence quark mass were generated by inserting the fit's determined 
values for its free parameters into ( |6.19| ) and ( |6.20| ). Recall that we use the pion- 



decay-constant normalization in which f.^ ^ 92.4 MeV. 

The diamonds display the result of a set of independent fits of the correlator at 



separate valence quark masses, using the same method as was used in Section pA . 
The result is a value for the pion mass and decay constant at each valence quark 
mass at which the correlator was calculated. The valence-quark-mass values at which 
a filled diamond appears correspond to the set of values used in the full pqChPT fit; 
that is, those within the valence-quark-mass cutoff. Open diamonds correspond to 
valence-quark-mass values beyond the cutoff. We stress that the full pqChPT fit is 
not a fit to the filled diamonds. Rather, the fit curve and points are related only in 
that they are derived from the same correlator data. Their agreement, or lack thereof, 
demonstrates the correlator data's tendency to match the predictions of pqChPT. 

Note that while two plots are used to present the results of the pqChPT fit, they 
are both the product of a single fit, encompassing the valence-quark-mass dependence 
of both the pion mass and decay constant, as well as the correlator's time dependence. 

Error bars were determined via a jackknife analysis of the statistical error. The fit 
curve's one-sigma range appears in the plots bound by dotted lines. This range was 
determined by individual jackknife analyses of the curve's value at each point along 
the horizontal axis. Uncertainty in the lattice spacing and kaon quark-mass threshold 
are not taken into account. Instead, the plots' axes have simply been rescaled using 
their central value, so as not to obscure statistical error in the quantities of interest. 

A diagonal correlation matrix was used in the fit, as the full matrix proved to not 
be positive definite in all cases. 
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Because the statistical error in Figure p.52| is too small to be visible, we have 



included, inset in the plot, close-ups at three values of the valence quark mass. These 
inset plots do not share a single scale. However, they do have a horizontal-to-vertical 
scale ratio which matches that of the main plot. 



The fit's determined values for its free parameters appear in Table |975 . 

9.6.3 Rm 



As is clear from Figure |D.52| , the dependence of the squared pion mass on the 
valence quark mass is very nearly linear. In such plots this strong linearity obscures 
the presence of higher-order effects. Thus, in order to accentuate non-linear terms. 



we include plots of a ratio first suggested in |^ 



_ myM^ms) 

where MT^^{mQ) represents the mass of the local pion containing valence quarks of 
mass iriQ. 

The ratio Rm is designed to allow a visual assessment of the strength and nature 
of the squared pion mass's non-linearity. Were the dependence of M^^ on my linear, 
with a massless pion in the chiral limit, a plot of Rm would be flat at Rm = 1- If the 
dependence were quadratic with no constant term, and did not include any higher- 
order or non-polynomial terms, an Rm plot would be linear. In both small 
non-zero pion mass at the chiral limit introduces a sharp downturn in Rm for small 
my- For larger my the behavior described above is unaffected. 



Taking the points and curve from Figure D.52 and transforming them based on 



the definition of Rm results in Figure D.54 
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In order not to cloud the statistical error of the fit results, a jackknife analysis 
of the quantity Rm has not been performed. As such, correlation between M^^ at a 
given valence quark mass and at the dynamical quark mass is not accounted for, nor 
is uncertainty in M^^ at the dynamical quark mass. Instead, the plot has simply been 
transformed using the central value for M^^{ms) as determined in Section |9.4] . We 



feel that this gives a more appropriate representation of the statistical error present 
in the pqChPT fit. The goal of presenting an Rm plot is not to accurately determine 
Rm and its uncertainty, but rather to compare the non-linearity of the data and 
the resulting fit curve. If we were to account for correlation with and uncertainty in 
[ms), it would obfuscate the points' error bars and the fit curve's one-sigma range. 
Furthermore, because the independent points and the fit curve do not necessarily agree 
on the value of M^^^{ms), a jackknife analysis would generate a misleading relative 
shift between their values for Rm- 

Note that we use Rm for plotting purposes only. For the pqChPT fit of the 
correlator data, the expressions for pion mass and decay constant, ( |6.19| ) and ( |6.2CI| ), 



are used directly. Additionally, we do not use the simplification of Rm suggested in 



9.6.4 Results 



Figures p.55| through |D]67| display the results of each ensemble's pqChPT fit. For 
each ensemble, three plots are presented: an Rm plot, a pion-mass plot, and a pion- 



decay-constant plot. The creation and presentation of these plots mirrors Figures D.52 



through p.54| , as discussed in Sections |9.6.2| and |9.6.3| . For completeness, the plots 
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ao 


z 


/ 


«5 


2^8 - "5 


A hyp 


0.025 


- 


9.35(17) 


0.05237(35) 


0.240(33) 


0.275(17) 


B hyp 


0.025 


- 


11.95(74) 


0.0471(11) 


0.857(58) 


0.211(53) 


C hyp 


0.025 


0.0979(77) 


149.1(85) 


0.00693(49) 


5.02(14) 


-0.28(11) 


A 


0.0125 


- 


9.60(15) 


0.07778(50) 


2.065(20) 


0.360(17) 


B 


0.0125 


- 


10.17(32) 


0.0758(11) 


2.187(49) 


0.356(41) 


C 


0.0125 


0.0411(41) 


43.7(20) 


0.0317(14) 


4.82(25) 


0.55(13) 


W hyp 


0.06 


- 


2.867(15) 


0.12036(26) 


-0.298(15) 


0.2472(66) 


X hyp 


0.06 


- 


3.104(19) 


0.11481(29) 


-0.066(21) 


0.2749(82) 


Y hyp 


0.06 


- 


3.281(23) 


0.11086(31) 


0.111(16) 


0.3004(90) 


Z hyp 


0.06 


- 


3.778(29) 


0.10221(32) 


0.491(17) 


0.3271(89) 


W 


0.055 




0.7649(20) 


0.23635(26) 


-0.527(26) 


0.278(14) 


X 


0.055 




0.8069(23) 


0.23049(27) 


-0.351(27) 


0.304(13) 


Y 


0.055 




0.8430(26) 


0.22548(28) 


-0.089(23) 


0.348(11) 


Z 


0.055 




0.9176(29) 


0.21596(28) 


0.254(25) 


0.4385(99) 


Q hyp 


0.025 




9.05(21) 


0.05250(46) 


-0.014(53) 


0.223(30) 



Table 9.5: Results from the pqChPT correlator fits and corresponding valence-quark-mass 
cutoff values. 



for ensemble A hyp are repeated. The values determined for the fits' free parameters 
are compiled in Table P3. 



Comparison of the pqChPT-predicted curve to the individual points in each en- 
semble's Rm plot demonstrates that the correlator data is systematically missing 
the predictions of pqChPT. This is most evident in the results for ensembles W hyp 
through Z hyp, and their thin-link counterparts. It is telling that these are the same 
ensembles which have the coarsest lattice spacing, and thus the largest expected 
flavor-symmetry breaking. It is likely that a majority of the data's systematic devia- 
tion from the predictions of pqChPT is due to our failure to account for the effects of 
the staggered formulation's inherent flavor symmetry breaking. As will be discussed 
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in Section 11.3 , any future work in this area will require a more robust handling of 
flavor-symmetry-breaking effects. 

Our results for the GL coefficient combination 2a^ — are relatively stable be- 
tween ensembles. Comparing ensembles A hyp, A, B hyp, and W hyp demonstrates 
that, while the corresponding systematic effects are strong, they are not beyond con- 
trol. A detailed analysis of the systematic error present in our calculation of 2a8 — 



is presented in Section [10.1.1 



Results for 05, which controls the polynomial NLO term in the pion decay con- 
stant, are much less consistent. This is interesting when we note that the data 
generally follow pqChPT's predictions for the form of the pion decay constant better 
than its predictions for the form of the pion mass. 

The determined values for 2^8 — from ensembles W hyp through Z hyp show a 
definite trend, demonstrating an apparent dependence on the dynamical quark mass. 
This trend can be explained if we recall the m^-dependent term which was dropped 
between the true predictions of pqChPT ( |6.17| ) and the form used in our single- 



dynamical-quark- mass fits ( |6.1!j| ). This term is accounted for in our simultaneous fit 



over these four ensembles, which is presented in Section 9.7 



9.6.5 Polynomial Fits 

For comparison we include the Rm plot for both a quadratic and a cubic polyno- 
mial fit of the squared pion mass of ensemble A hyp. These can be found in Figures 
D.68| and p.69| . The quadratic fit is taken directly from Section |9.5| , while the cubic 



fit uses the same methods as Section |9.5| , replacing the fit's form for the squared pion 
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ao ai a2 03 


A hyp 


0.001576(82) 3.815(28) -8.12(81) 78.0(81) 



Table 9.6: Results for a subset of the fit parameters from the cubic fit of the pion mass's 
valence-quark-mass dependence. 



mass (19. 61) with: 



M = ao + aiTTLv + a2my + a^my (9.9) 



The resulting values for the fit parameters of the cubic fit are given in Table ^ 
The agreement between the independent mass points and the fit curve, even in the 
quadratic case, is strikingly good, especially when compared to the results of the 



corresponding pqChPT fit. Figure D.54 



9.6.6 Finite Volume 

The fit results of ensembles C hyp and C are given last, as their analyses re- 
quired special attention. The small volume of these ensembles significantly affects 
their correlator data, a fact which is made most clear by the large non-zero value 
obtained for the constant term in the calculation of their kaon quark-mass thresh- 
olds. The reader is directed to Table |9.4| and Figure |D.47|. These finite- volume effects 



overwhelm the forms predicted by infinite-volume pqChPT, such that attempts to fit 
these ensembles' correlator data to pqChPT's forms either fail completely or generate 
nonsensical results. 

In order to produce an even moderately reasonable fit, we added a constant term 
to pqChPT's predictions for the form of the squared pion mass. Thus, in these 
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fits, the standard pqChPT form ( |6.19| ) was replaced with: 

Ml^ = ao + zmy(47r/)2|l + — (2my - ms) In zmy + Jf^i^v - mg) 

+ zmv{2as-a5)\ (9.10) 



where ao is an additional fit parameter. Comparing Table |9.4] and |9.5| , it is interesting 
to note the similarity between the values obtained for the constant term oq in the 
quadratic and pqChPT-predicted forms. 



Even after the addition of a constant term, the results of the fits. Figures p.7C 



and |D.71|, remain questionable, especially in the case of ensemble C hyp. Luckily, our 



study of these ensembles has no purpose other than to elucidate the effects of finite 
volume. 

9.6.7 Dependence on Amy 

In order to study the dependence of our results for 2as — as on the choice of 
valence-quark- mass cutoff Amy, the pqChPT fits were repeated using a range of cutoff 
values. The results of this investigation are presented in Figures p.72| through p. 76 



Figure p.72| displays the A^y dependence of all four parameters of the fit of ensemble 



A hyp, while Figures p.73| through p. 7^ display only the dependence of the parameter 



2as — as for all ensembles. 

These plots clearly demonstrate that the cutoff choice is a significant source of 
systematic error. The strong dependence of 2a8 — as on A^y is a direct result of 
the failure of the theoretical forms to closely match the data. As the mass cutoff is 
increased, and additional valence-quark-mass values are added to the fit range, the 
additional correlator data consistently fails to match the values predicted for it by 
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previous lower-cutoff fits. As such, each addition of correlator data significantly alters 
the fit results. 

We expect this sort of behavior for large cutoffs which fall beyond the range 
of pqChPT. However, below some threshold value for A^^, we expect the value of 
2a8 — as to level off. This threshold would indicate the outer limit of pqChPT's 
domain. Yet, in our data we see no clear plateau. Instead, we observe a strong 
dependence on K^y down to the smallest valence quark masses studied. 

Two possible reasons for this behavior are readily available. First, it is possible 
that the valence quark masses under study are too large for pqChPT to generate 
accurate predictions. That is, our entire study lies beyond the threshold. However, 
as the mass of our local pion reaches values well below the physical kaon mass, we feel 
that this possibility is unlikely. Second, the behavior of our correlator data may be 
significantly skewed by fiavor-symmetry-breaking effects, rendering the predictions of 
continuum pqChPT inaccurate. 

While the choice of valence-quark-mass cutoff is clearly a significant source of 
systematic error, the variation of 2a^ — is not so great as to render our results 



meaningless. In Section |10.1.1| we incorporate this quantitative analysis of this source 
of systematic error into an estimate for the full systematic uncertainty of our result. 

While the utility of an analysis of the A^^ dependence of the results of ensembles 
C hyp and C is not clear, it has been included for completeness. Because, for these 
ensembles, we are using a pion mass form with three free parameters, the smallest 
cutoff which generates sensible results is one which leaves us with three valence-quark- 
mass values within the cutoff. As such, the range of studied for these ensembles 
does not reach as low as for other ensembles. 
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/ 


2^8 - «5 


"5 


2^6 ~ Q^4 


0^4 


hyp 
thin 


2.500(16) 
0.6970(26) 


0.12993(41) 
0.24780(39) 


0.3860(63) 
0.2292(88) 


-0.119(19) 
-0.093(29) 


-0.2703(74) 
0.006(14) 


-1.446(33) 
-3.207(59) 



Table 9.7: Results from the simultaneous pqChPT fits of the correlators of ensemble sets 
hyp, which includes ensembles W hyp through Z hyp, and thin, which includes ensembles W 
through Z. 



9.7 2^6 ~ 



Four of our ensembles, ensembles W through Z, are constructed to have the same 
lattice spacing and volume, such that the only variation between them is their dy- 
namical quark mass. Through these ensembles we are granted the opportunity to fit 
lattice data to the predictions of pqChPT over a block of pqQCD's two-dimensional 
quark-mass plane. 

We fit the correlator C^^^-^t over a range of time separations t, valence quark masses 
my, and dynamical quark masses ms to the form predicted by pqChPT: ( 6.16 ), ( |6.17| ), 



and ( |6.18| ). Because we are working over a range of dynamical quark masses, we do 



not need to drop the unknown ms dependence from these forms, as was done in 



Section |9]^. The result is values for the fit's free parameters, which include the GL 
coefficient combinations la^ — 0:5, — 04, as, and 04. 

The fit was carried out twice, once for the hypercubic-blocked ensemble set, and 
once for the corresponding thin-link set. The resulting fit parameter values can be 



found in Table |9.71 , with the ensemble set hyp containing ensembles W hyp through 



Z hyp and the ensemble set thin containing ensembles W through Z. The correspond- 



ing fit curves are displayed in Figures p.77| through p.80| . These curves are generated 
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using the fits' resulting parameter values in equations ( |6.17|) and ( |6.18|) and construct 



ing the Rm curves as described in Section p.6.3| . Each plot includes four cross sections 



through the quark-mass plane, each along a different line of constant dynamical quark 
mass. Note that while eight fit curves are shown for each ensemble set, four in each 
plot, they together represent the results of a single fit. 

The plots' diamonds represent the results of a set of independent fits of the cor- 
relator data at separate valence- and dynamical-quark-mass values, using the same 
method as was used in Section From this we obtain a value for the pion mass 
and decay constant at each valence and dynamical quark mass at which the correlator 
was calculated. The quark-mass values at which a filled diamond appears correspond 
to the set of values used in the full pqChPT fit. Open diamonds correspond to 
quark-mass values beyond the valence- quark- mass cutoff. The full pqChPT fit of an 
ensemble set is not a fit to the filled diamonds. Rather, the fit curve and diamonds are 
related only in that they are derived from the same correlator data. Their agreement, 
or lack thereof, demonstrates the correlator's tendency to match the predictions of 
pqChPT. 

Error bars were determined via a jackknife analysis of the statistical error. The fit 
curves' one-sigma range appears in the plots bound by dotted lines. This range was 
determined by individual jackknife analyses at each point along the horizontal axis. 
Uncertainty in the lattice spacing and kaon quark-mass threshold are not taken into 
account. Instead, the plots' axes have simply been rescaled using the average central 
value between the ensembles. 

A diagonal correlation matrix was used in the fit, as the large number of correlator 
values generated a matrix which was far from positive definite. 
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For illustrative purposes we have compiled the results of the independent pqChPT 



fits of ensembles W hyp through Z hyp and W through Z from Section into plots 
whose formats mirror Figures p.77| through p.80| . These compiled plots of the inde- 
pendent fits can be found in Figures U.81 through U.84 . Comparison between these 
plots makes clear the difference in results between the independent and simultaneous 
correlator fits. 

We also present the results of the simultaneous fits using a cross section through 
the quark-mass plane along the unquenched line, mq = my = ms- The plots are 



found in Figures p.85| and |D.86| . Other than the choice of cross section, these plots 
were generated in the same fashion as described above. Because the definition of Rm 
is meaningful only along lines of constant dynamical quark mass, we introduce a new 
ratio: 



For our plots we use a reference quark mass of itir = 0.025. 



(9.11) 



From all the plots presented, we can clearly see that the correlator data system- 
atically misses the predictions of pqChPT. As discussed in Section |9.6.4| , this is most 
likely due to strong flavor-symmetry-breaking effects, a consequence of these ensem- 
bles' coarse lattice spacing. Because of the computational expense involved in the 
generation of a set of four reasonably long partially quenched Markov chains, we were 
forced to use small lattice extents. Thus, a very coarse lattice spacing was required 
in order to maintain a reasonable lattice volume. The spacing is coarser than what is 
generally deemed acceptable by the community. As such, this calculation of 2aQ — 
can only be taken as a preliminary study. Yet, as ours is the first attempt at such a 
calculation, such a preliminary study is not without value. 
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The preliminary nature of this study is further emphasized by the fact that only a 
single ensemble set is available to us. As such, we do not have the ability to generate 
estimates for the magnitude of the various systematic errors we know are present in 
our result. 



9.7.1 Dependence on A 



my 



In order to study the dependence of our results on the choice of valence-quark- 
mass cutoff Am^,, we repeated the fits using a range of cutoff values. The results 
of this investigation are shown in Figures p.87| and p.88|. Just as was seen in the 



independent fits of Section p.6| , there is a strong dependence on Kmy down to small 
quark mass. We would hope that in a more complete study using a smaller lattice 
spacing, this dependence would disappear for small quark mass. 
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CHAPTER 10 



RESULTS 



The primary result of our study is a value for the Gasser-Leutwyler coefficient 
combination 2Ls — L5 (|10.3|) , along with its corresponding value for the light-quark- 
mass ratio rriu/md ( |10.9| ). The secondary results include a value for the the GL 
coefficient L5 with very large systematic errors ( |10.11| ) and values for L4 and Lq with 
large and unestimated systematic errors, ( |1U.12| ) and ( |1(J.13D . 



A subset of these results were presented in |jT|. Their account here is signifi- 
cantly more comprehensive and incorporates several improvements in our analysis 
techniques. 

10.1 Primary Results 

We take the central value of our quoted result for the GL coefficient combination 
2as — «5 from the correlator data of our primary ensemble, ensemble A hyp. This 
produces the value 2as — 0^5 = 0.275 ± 0.017, where the given uncertainty accounts 
only for statistical error. 

10.1.1 Systematic Error 

An important aspect of our study is the ability to make a quantitative estimate 
of our systematic error. Past theoretical estimates for the GL coefficient combination 
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2as — a5 have suffered from an inability to quantify the systematic error resulting from 
the approximations they require. Our first-principles approach, however, allows for 
such an investigation. All of the systematics which separate our calculation's result 
from the true value of 2a8 — are known and open to investigation. 

To determine our systematic error, we note the variation in 2as — ck^ due to four 
changes in our calculation: reducing the lattice volume, hypercubic blocking, doubling 
the lattice spacing, and shifting the valence-quark-mass cutoff. While it is likely that 
several of the variations are correlated — for example, both the removal of hypercubic 
blocking and doubling the lattice spacing increase a single systematic-error source: 
flavor and Lorentz symmetry breaking — in order to generate a generous estimate of 
our error, we will add the variations in quadrature, as if uncorrelated. 

Our study of ensemble B hyp grants us insight into the effects of finite volume on 
our result. Ensemble B hyp produces values for the fit parameters similar to those 
of ensemble A hyp, indicating that finite-volume effects in our primary ensemble are 
well under control. We estimate the uncertainty due to finite volume using the shift 
in value of 2as — between the volumes: ±0.064. 

Hypercubic blocking has been shown to reduce the flavor symmetry breaking in- 
herent in staggered calculations. In order to estimate the strength of flavor-symmetry- 
breaking effects on our result, we compare our result to the value obtained from the 
thin-link version of our primary ensemble, ensemble A. We take the resulting shift in 
value as our estimate of flavor-symmetry-breaking error: ±0.085. 

The lattice spacing has a direct impact on the accuracy of our lattice calculations, 
controlling the strength of unwanted Lorentz- and fiavor-symmetry-breaking terms in 
our action. To estimate the strength of these finite-lattice-spacing effects, we compare 
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our value for 2a^ — to that obtained from an ensemble with approximately double 
the lattice spacing, ensemble W hyp. This shift gives us a qualitative estimate for the 
corresponding systematic error: ±0.028. 



The investigation of Section p.6.7| demonstrated that the choice of valence-quark- 
mass cutoff is a significant source of systematic uncertainty. While this uncertainty 
is likely related to flavor symmetry breaking, we will add in quadrature the variation 
due to changing our valence-quark-mass cutoff as if it were an independent error 



source. It was determined in Section |9.6.7| that shifting the valence-quark-mass cutoff 



^mvl^QK by ±0.15 results in a variation of 2as — 05 equal to ±0.13. We will use 
this variation as our estimate of the systematic error due to our somewhat arbitrary 
choice of valence-quark-mass cutoff. 

Presenting all error sources together, our result becomes: 

2as-a^ = 0.275 ± 0.017 ± 0.064 ± 0.085 ± 0.028 ± 0.13 (10.1) 

Adding all sources in quadrature gives our final result: 

2^8-05 = 0.28 ±0.17 (10.2) 

In terms of the standard GL coefficient normalization, this corresponds to: 

2^8 - L5 = (0.22 ±0.14) X 10"^ (10.3) 



Note that this lies outside the range in which a massless up quark is allowed ( |3.42| ) 



The compiled results for 2^8 ~ ct5 from all ensembles, as well as both our final 



result and the range which allows a massless up quark, can be found in Figure |10.1 
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Figure 10.1: The compiled results for 2as — 05. The bold open diamond corresponds to 
our final result and our quoted statistical and systematic error. All other error bars are 
statistical. The gray diamonds correspond to the result from ensemble A hyp after shifting 
the valence-quark-mass cutoff /f^Qx by ±0.015. The filled and open circles correspond 
to the results of the simultaneous fits of ensembles W hyp through Z hyp and W through Z 
respectively. The burst corresponds to the range of values which allow a massless up quark 
( 3.44 ). All other points represent the values obtained for 2a^ — 05 from their corresponding 
ensemble. 
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10.1.2 Light-Quark-Mass Ratio 



As discussed in Section p^ , the GL coefficients of our partially quenched calcu- 
lation are the same coefficients that appear in the physical chiral Lagrangian. As 
such, we can use our partially quenched results for 2L^ — L5, along with NLO ChPT 
as presented in Chapter to determine the light-quark-mass ratio of the Standard 
Model. 

Using our result ( |10.3|) in ( |3.35D generates a value for the NLO correction Ajv/: 



0.0919 ± 0.029 ± 0.0084 (10.4) 



where the ffist uncertainty is due to error in our result and the second is based on 
an assumption that the unaccounted for NNLO corrections are on the order of A|^. 
Error in other inputs to the calculation of Am are overwhelmed by the uncertainties 
given. Using Dashen's Theorem to account for the QED contributions to the physical 
meson masses, instead of ( p. 291 ), does not affect Am at this level of precision. 



Using our result for 2L^ — L5 ( |10.3| ) in ( p.40| ) allows us to generate a value for the 
light-quark-mass ratio: 



Tfl 

— = 0.408 ± 0.027 ± 0.008 ± 0.021 (10.5) 



where the ffist uncertainty is due to error in our calculated result, the second comes 
from an assumption that the NNLO corrections to Am are on the order of A\.j, and 
the third is due to uncertainty in A^;. 
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Summarizing the repercussions of our calculated value for 2L^ — L^: 



Am = 0.092 ± 0.030 
2L8 - = (0.22 ± 0.14) X 10"^ 
(-0.23 ± 0.09) X 10"^ Ls = (0.36 ± 0.24) x 10^^ 

TTl 

— = 0.408 ±0.035 



(10.6) 
(10.7) 
(10.8) 
(10.9) 



where the error given in (|10.4|) and (|10.5|) has been added in quadrature. In order 
to determine L7 and Lg, we have used the value for L5 determined from the physical 
meson-decay-constant ratio ( |3.51| ) and the value for I2L7 + GLg — L5 determined 
from the physical deviation from the Gell-Mann-Okubo relation ( p.55| ), adding in 
quadrature the experimental uncertainty in those values to the error in our calculation. 
Note that we do not use our calculated value for L5 as presented in Section |10.2.1. 



These numbers can be compared to the generally accepted set of values, ( p.91| ) 
through (|3.94| ), which are obtained using various model-dependent assumptions and 
data from beyond the light-meson sector. 

The light-quark-mass ratio which results from using Dashen's theorem, instead of 
( p.29| ), is rriu/md = 0.482 ± 0.026, where there is now no uncertainty due to Ag. 



Our result for the light-quark-mass ratio (|10.9| ) can be compared with the pre- 



diction generated by the generally accepted value for Am (|3.91| ): rriu/md = 0.608 ± 
0.056. The literature provides additional predictions. Leutwyler gives a ratio of 



^u/f^d = 0.553 ± 0.043, while Amoros et al. [|T9| give a ratio of rriu/md = 0.46 ± 0.09 
( p.99| ). Our result is similar to these values, yet somewhat smaller in each case. Note 
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that our error bars are smaller than all quoted cases. In addition, we are of the opin- 
ion that our error bars are the first that can be well trusted, as all others are attempts 
to account for uncontrolled systematic error due to theoretical assumptions. 

The Review of Particle Physics |^8[ quotes a very broad range for the light-quark- 



mass ratio: 0.2 < mu/ma < 0.7. We fall well within this range. 

Our first-principles calculation of the low-energy constants of ChPT demonstrates 
that the relevant coefficients are too low to allow for the scenario in which the up quark 
is massless and strong NLO terms emulate a non-zero mass. As the massless-quark 
solution is seen here to be unlikely, the strong CP problem remains unanswered. 

10.2 Secondary Results 

The secondary results of our study include values for the GL coefficients L5, L4, 
and Lq. 

10.2.1 L5 

Our results for the GL coefficient L5 are extremely vulnerable to systematic error 
and vary significantly across the ensembles studied. Estimating the systematic error 
using the same method as Section |10.1.1| results in: 



as = 0.240 ± 0.033 ± 0.62 ± 1.8 ± 0.54 ± 0.008 (10.10) 

where the listed uncertainties are due to statistical error, reducing the lattice volume, 
hypercubic blocking, doubling the lattice spacing, and shifting the valence-quark-mass 
cutoff. While it is not clear with such large variation that adding in quadrature is 
justified, doing so results in = 0.2 ± 2.0. In terms of the standard GL coefficient 
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normalization, this corresponds to: 

L5 = (0.2 ± 1.6) X 10^^ (10.11) 

This value falls within the range obtained for L5 from the meson-decay-constant 
ratio ( p3l| ). 

10.2.2 L4 and Lq 

Due to the coarse lattice spacing of ensembles W through Z, our results for the 
GL coefficient combinations 2aQ — and 04 are guaranteed to be contaminated by 
strong systematic error. However, with only a single ensemble set to work with, 
we are unable to make any qualitative estimates of the error. Thus we present our 
values for these coefficients only as preliminary results, using the values from the 
hypercubic-blocked ensemble set and quoting only their statistical error bars: 

L4 = (-1.145 ± 0.026^"*^*)) X 10"^ (10.12) 
Lq = (-0.680 ± O.OIS^"*^*)) X 10-3 (10.13) 

These values both miss their generally accepted ranges, L4 = (—0.3 ± 0.5) x 10^'^ 
and Lq = (—0.2 ± 0.3) x 10"^, which are determined using large-A^^c considerations 

10.3 Quenching Effects 

Analyzing the results produced by ensemble Q hyp, we see that quenching had 
a smaller systematic effect on 2as — 05 and 05 than any other source of systematic 
error studied. Clearly, at this level of precision, quenched and partially quenched 
ensembles do not generate distinctly different valence-quark-mass dependencies for 
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the pion mass and decay constant, despite the differing predictions of quenched and 
partially quenched ChPT. This indicates that, while the Nf dependence of ChPT 
is not explicit and thus the GL coefficients are unknown functions of the number of 
dynamical quark flavors, that functional dependence is very slight. 
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CHAPTER 11 



SUMMARY AND OUTLOOK 



The low-energy constants of the chiral Lagrangian, the Gasser-Leutwyler coeffi- 
cients, are a critical element in the understanding of low-energy QCD. Yet current 
theoretical estimates and experimental measurements of the GL coefficients, even 
those unaffected by the Kaplan-Manohar ambiguity, have errors from 10% to 160% 
p3|. For many of the coefficients, their current error bars remain as large as they 



were at the ffist instance of their calculation [|T^. Given LQCD's capacity to calculate 
these coefficients directly, with no uncontrolled approximations, it is clear that work 
in this area is warranted. This may in fact prove to be one of those rare, yet increas- 
ingly common, situations in which lattice techniques have some chance of providing 
the greater community with the most trusted predictions available. 

Our study definitively calculates a single combination of the GL coefficients, de- 
termining a value for 2L^ — L5 which rules out the massless-up-quark solution to the 
strong CP problem. The culmination of our study is a value for the light- quark-mass 
ratio: 

Tfl 

— = 0.408 ± 0.035 (11.1) 
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This work is far from the closing word in the lattice study of Gasser-Leutwyler 
coefficients. In truth it is just the first step in what will likely be a long-term and 
comprehensive study of the coefficients by the lattice community. 

From a short-term and more pragmatic perspective, there are several aspects of 
our study ripe for improvement. In addition, new theoretical work is on the horizon 
which should dramatically improve lattice calculations of the GL coefficients. 

11.1 Improved Systematics 

Clearly, future calculations would do well to improve on the systematics of our 
study through superior ensembles, utilizing either more sophisticated actions or sim- 
ply more sophisticated computers. 

Despite the large uncertainty in the QED contribution, the greatest source of error 
in the light-quark-mass ratio (|10.5|) remains the systematic error in our calculation 
of 2L^ — L5. Improved systematics thus have the potential to reduce the error bars 
by up to a factor of root two. 

Additionally, a reduction in systematic error may prove to bring the fluctuations 
in L5 under control, leading to a result in which we could have reasonable confldence. 
As L5 is not subject to the KM ambiguity, this would allow for a comparison between 
lattice and experimental results. 

In the case of the coefficients L4 and Lg, our study was not broad enough to 
estimate the magnitude of our error. Producing results with reasonable and well- 
estimated systematic error would require the availability of a number of quality en- 
semble sets, with each set including ensembles across a range of dynamical quark 
masses. 
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11.1.1 High-Quality Publicly Available Ensembles 

When our study began, high-quality ensembles at Nf = 3 were not publicly avail- 
able. We had no choice but to use our relatively limited computer resources to 
generate the requisite ensembles. As a consequence the resulting ensembles, with 
their mediocre lattice extent and unimproved action, are far from cutting edge. 

Today however, thanks to the MILC collaboration |7T| in conjunction with the 



Gauge Connection [Q, a set of ten 20^ x 64, Nf = 3 ensembles are available for public 



use. These ensembles span a range of dynamical quark masses and were generated 



using an improved gauge action and the highly improved a -tad staggered action [[73 
The application of our analysis to these ensembles would not only produce reduced 
systematic error in the calculation of 2Lg, — L5, but as the ensembles have matched 
lattice spacings, would also allow for an accurate calculation of L4 and Lq. 

Such a spirit of sharing is extremely valuable to those of us who do not lead the 
community in computational resources. Beyond that, the community itself profits, as 
the number of individuals capable of significant and impacting research is dramatically 
increased. 

11.2 Dynamical Hypercubic Blocking 

In our study we reduced the systematic error due to flavor symmetry breaking 
by hypercubic blocking after the generation of our ensembles. In effect we used 
hypercubic-blocked valence quarks in conjunction with thin-link dynamical quarks. 
While this technique has merit, clearly a more consistent approach would be to 
use hypercubic-blocked dynamical quarks during ensemble creation. This approach 
should also further decrease fiavor-symmetry-breaking error. 
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Following this logic, Flemming [f^] has plans for a direct continuation of our work 
via the creation of ensembles using practical algorithms for dynamical hypercubic- 
blocked fermions, which are only now becoming available [KO, |51|, |5^ . 



11.3 Staggered Chiral Perturbation Theory 

A recent and very exciting theoretical development germane to the connection 
between LQCD and ChPT is staggered Chiral Perturbation Theory (sChPT) ^ . 

In standard ChPT only a single structure breaks flavor symmetry: the quark 
mass matrix. In the context of staggered fermions, however, the theory contains 
additional flavor-breaking structures which arise at finite lattice spacing. In sChPT 
these new flavor-breaking elements are accounted for through the introduction of 
additional terms to the Lagrangian at each order. The chiral Lagrangian becomes an 
expansion in three parameters instead of two: meson momentum, meson mass, and 
lattice spacing. 

A complication arises in sChPT when work is done at Nf ^ 4. At the core of 
the staggered formulation are four flavors. Thus, even when a fractional power of 
the fermionic determinant is used to set Nf ^ 4, as discussed in Section 5.3.9 , the 



flavor symmetry breaking due to finite lattice spacing retains its four-flavor structure. 
As a consequence, sChPT can only be constructed for theories which include some 
multiple of four flavors. In order to apply the results of sChPT to theories with some 
other number of flavors, such as Nf = 3, the strength of meson-loop graphs must be 
adjusted by hand. 

As discussed in Section ^.3.6| , flavor symmetry breaking in staggered fermions 
leads to a non-degeneracy of the sixteen light mesons, splitting them into five levels. 
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The sChPT expression for the local pion mass conforms to this scenario, containing 
additional terms which are the result of meson loops failing to cancel exactly, as they 
would in the degenerate continuum case. As such, Bernard W^ demonstrates that 



accounting for a majority of the new finite-lattice-spacing terms does not require the 
calculation of a full set of fresh low-energy constants. Instead, it requires only a 
calculation of the four meson-mass splittings. 

However, two of these new terms, those which allow mesons to switch flavor con- 
tent via insertions into a propagator, have coefficients which can not be directly 
measured. Thus, they must be left as free parameters in a fit. Preliminary attempts 
at such fits of the local pion mass have proved unstable, with the two parame- 



ters running to unnaturally large and opposite values. Introducing priors to the fit 
could stabilize results. Also, as the same two coefficients appear in all sChPT ex- 
pressions, fitting multiple quantities simultaneously, such as the pion mass and decay 
constant, may bind their values. It is worth noting that these coefficients can not be 
determined once and for all, unlike the GL coefficients, as they are a function of the 
action. Different improved actions have different flavor symmetry breaking, and thus 
the coefficients will take on different values. 

Because the lattice-spacing dependence of the staggered chiral Lagrangian is known 
and explicit up to whatever order we chose, and because in the continuum limit stag- 
gered quarks are equivalent to continuum quarks, a calculation of the low-energy 
constants of sChPT corresponds exactly to a calculation of physical ChPT's GL co- 
efficients. 



Preliminary results for sChPT are promising |^. They demonstrate that in 
cases where the local pion mass's dependence on the quark mass does not follow the 
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form predicted by continuum ChPT, a phenomenon clearly evident in our data, the 
dependence matches closely the predictions of sChPT. As such, use of sChPT should 
result in a drastic reduction of the systematic errors in lattice calculations of the GL 
coefficients. 

In fact the advent of sChPT has an impact beyond simply the accurate measure- 
ment of GL coefficients. sChPT provides, for any quantities to which it is relevant, 
the appropriate simultaneous extrapolation from numerically-favorable quark masses 
and finite lattice spacing to physical quark masses and the continuum limit. 

11.4 Quantum Electrodynamic Corrections 

A significant percentage of the uncertainty in the light-quark-mass ratio is due 
to uncertainty in the magnitude of the QED contribution to the pion mass. While 
not addressed in this study, lattice techniques exist which allow for the calculation of 
this contribution [|77|, |T8[. In a situation where attempts to significantly reduce the 
systematic error in 2L8 — L5 are successful, an accurate lattice calculation of A^; may 
prove valuable. 
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APPENDIX A 



NOTATION 



So as to leave the main text uncluttered, we list several notational conventions 
here. 

Throughout, we use the pion decay constant normalization /jr — 92.4 MeV. This 



When the standard normalization for the GL coefficients is being used, they are 
denoted by Lj. While this is the normalization which generally appears in a chiral 
Lagrangian, most NLO expressions for observable quantities are made cleaner through 
the use of a second normalization which we denote by a^. The normalizations are 
related by the expression: 



Traces over color indices are denoted by tr, while traces over other indices, gener- 
ally flavor indices, are denoted by Tr. 

Integration over all space is denoted by: 



differs from the other common normalization by a factor of root two, V^/tt — 130.7 MeV. 



(A.l) 




(A.2) 



while integration over all momenta is denoted by: 




(A.3) 
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When denoting dimensionful quark masses, we use a lower-case subscript such 
as rriq. When denoting a dimensionless lattice quark mass, we use a corresponding 
upper-case subscript such as mg. For all quantities other than quark mass, we use a 
check over the variable to denote its unitless counterpart. For example in the case of 
the local pion mass, M^g = o-^ttb- 

The generators of the various SU{N) Lie algebras are represented by A" for color 
transformations and r" for flavor transformations. They are traceless Hermitian 
N matrices, which are normalized according to: 

TrfrV^ = ^5"^ tr[A"A''] = ^5"^ (A.4) 

and satisfy the commutation relations: 

[r«, T*"] = [A", A*"] = ig^'^'X" (A.5) 

where /"^"^ and g'"'"^ are the appropriate structure constants of the algebras, which 
are completely antisymmetric and real. A" should not be confused with the SU{N) 
Gell-Mann matrices, as A" = ^X.Q^yy_^^^^. 

The Lorentz tensor e/j^i^a/s is defined to be totally antisymmetric, with: 

61234 = 1 (A.6) 

The Euchdean-space Dirac matrices 7'* are defined to satisfy the anticommutation 
relation: 

{7^ 7"} = 2^^^ = 2(5'''^ (A.7) 

where g'^ is the fiat Euclidean-space metric. Unlike the Minkowski-space Dirac ma- 
trices, the Euclidean-space matrices are Hermitian: 

Y = 7^t (A.8) 
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The Dirac matrices can be used to construct the generators of a spinor representation 
of the Lorentz group. Given a Lorentz transformation: 



^ x'^^A.i'X (A.9) 

K = e-t-"'^(^"'')^^ e 0(4)Lorentz (A.IO) 

where J"'^ are the generators of Lorentz rotations on Lorentz 4- vectors and ujai3 is 
an antisymmetric tensor which parameterizes the transformation, a spinor transforms 
as: 

^ ^ il}' = S{K)il} (A.ll) 
^ ^ ^' = V55-i(A) (A.12) 
S{K) = e-5-«M^'".T''l G 0(4),pi„ (A.13) 

The Dirac matrices themselves transform correctly under the Lorentz group both as 
Lorentz 4- vectors and as spin space matrices: 

5(A)7^5^i(A) = K'^.Y (A. 14) 

Using the Dirac matrices, the matrix 7^ is defined as: 

75 = 7S'7'7^ (A. 15) 

which anticommutes with the other four Dirac matrices: 

{75,7/^} = 7Y = 1 (A.16) 
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We define the discretized directional derivative as: 



2a 



-fsinh ad^}\ip{x) 

d^ip{x) + 0{a^) (A.17) 



while we defined the second discretized directional derivative A^, as: 



Ai^ix) 



(fi{x + afi) — 2(f{x) + (f{x — afi) 



0? 



1 

a? 

= 4 (l + + ¥^^1 - 2 + 1 - + \a'd''^^{x) + O(a^) 

^ dl^{x) + 0{a') (A.18) 
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APPENDIX B 



STAGGERED IDENTITIES 



B.l Bilinear Sum 

For the sources of our bilinear correlators, we use a linear combination of all 
bilinears with a given distance binary four-vector T>s,j^ — T>. Such a combination has 
a simple form when expressed in terms of % and 



Jn+v,n;h = ^ Qh{lv1n ® ^Ti)Qh 
n n 

n 

= 4Tr[g;,7^]Tr[Q^] 

= XI XI [r^7x,] Tr [Fb] Xh+AXh+B 

A B 

= 16 X X ^A,'D5B,0Xh+AXh+B 



A B 

= l&Xh+vXh (B.l) 
where we have used the identity 

E(ra.,(rc)„„ = 4MAa (B.2) 

c 

or equivalently 

Y,^\VcTbTI = AV\TtVb (B.3) 
c 
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The linear combination leaves only a single contraction, that between x at the corner 
of the hypercube offset from the lowest corner by V and x fhe lowest corner. 
Adding the gauge links required by an interacting theory, we have: 



(B.4) 



n 



where lAA,B;h is defined after ( ^.621) . In the case of local bilinears, P = 0, the contrac- 
tion is contained completely by the lowest corner of the hypercube: 



^ Jn,n;h = 16x, 



hXh 



(B.5) 



n 



B.2 Transpose of the Staggered Interaction Matrix 

The transpose of the staggered fermion interaction matrix M^[U]^ arises as the 
antiquark propagator in our calculation of bilinear correlators. The adjoint matrix 
can be expressed in terms of M'^[U] in a manner similar to that used for the naive 
interaction matrix (|5.28| ): 



fi;n n,m— fl fi;n— fi n;m+ fi 



(B.6) 



This mirrors closely the result for naive fermions, with e„ filling the roll of 75. 
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B.3 Non-Local Bilinear Correlators 



In order to study the non-local staggered mesons, the correlators between non-local 
bilinears must be calculated. We now express non-local bilinear correlators explicitly 
in terms of quantities which are straightforward to calculate using lattice techniques. 
That is, the inverse staggered fermion interaction matrix applied to various fermion 



field vectors. A similar discussion limited to local bilinears is found in Section [5.8.1 . 

Our general bilinear correlator, using a wall sink to overlap only with zero mo- 
mentum states, is: 



9 

We replace the bilinear in our source with a linear combination of all bilinears with 
distance vector V = Vs jr = S + 

Y^(y^ Js,r;ay^,Jn+v,'R.;o) (B.8) 

g 7^ 
g4=t 

We replace our single-bilinear source with a wall of bilinears at the appropriate time 
slice: 

g fl TZ 

gi=t h4=0 



Using (|5.62| ), ( |5.63| ), and ( |B.4| ), we express the correlator in terms of x and x'- 



9 B 



E E Xg+BUB,B+V;gXg+B+VXh+vUvfl-hXhj (B.IO) 



_ , , a a,b b d d,c 

h a,,b,c,d 

hA=0 
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Integrating over the Grassmann variables using Wick contractions, we express the 
correlator in terms of the inverse of the interaction matrix: 



9 B 
g4=t 



X X] X9+BUB,B+V;gXg+B+VXh+vUvfi;hXh 

, J a a.h b d d.c c 

h a,b,c,d 

h4=0 



= (EEh 



9 B 

h a,b,c,d "'^ '''^ "''^ b,d 

h4=0 

9 B 

E E ^s.^^+^;5^^.o;/.M^[C/];^5,,M^[C/];^^+^,,+^) (B.ll) 

h a,b,c,d '^''^ "''^ ''''^ 

/l4=0 

In order to reduce the number of required applications of the inverse interaction 
matrix, we calculate the 2N(. field vectors X^'^^ and Y^'^^: 

XL^^ = ^ M^[C/]i (B.12) 



a,c 
ft 

ft4=0 



= E E ^''[^]-!/^+I^W^',0;/^ (B.13) 



/l4=0 



To calculate X^^^ we construct the field vector W^'^\ which equals one only at color 
c in the lowest corner of each hypercube on the time slice 714 = and zero elsewhere. 
The result of applying the inverse interaction matrix to W^'^^ is X^"^); 

xiC^ = {M^[U]-'W^% (B.14) 

a a 
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To calculate Y^'^^ we first apply a swap operator Sx> to W^'^^ and then the inverse 
interaction matrix. The resulting field vector is Y^'^^: 



yif) = {M'^[U]-^ST,W^'^)m (B.15) 

b 



The swap operator Sx> is similar to the bilinear operator (75 ® defined by (|5.75| ), 
but does not apply phase factors to the field vector. It swaps corners separated by 
the offset vector V within each hypercube, and applies the appropriate color matrix 
for that movement: 

('^»)n,m = (B.16) 

where h denotes the hypercube containing n, and A denotes its position within that 
hypercube: 



h = 2 



A = n-h (B.17) 



The unique random phase associated with each lattice site due to local gauge 
freedom washes out the position-off-diagonal terms in the product of X^'^'' and F*-^^: 



b a ^ b^d d,c ^ a,c 

h f 
h4=0 /4=0 



= E E M''[U]-]l^jyl(v,o;hM'[U]-X (B.18) 

■ ^ b.d d,c a.c 

d h 
/i4=0 

Thus, once X*^^) and y'-'^-' have been calculated, we can construct our correlator: 

cs,., = (E E E(-)"^^^'^^^" Y.^b,b.-,y;^%,^x%) 



c g B 
94=t 



c g B a,b 
94=t 



E E E(-)'^"^'""^'" E^15i^^^.^+^^^n%+i,) (B.19) 

a,b f, I 
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or more concisely 

cs,^;t = (E E E(-)"'"^'""^^" E^i?B(^;).+i^,.+B+.,i;%+^) 

c g B a,b "'^ b ' 

j:i-r^''^^'-'''Y.^lUsvY^%^B) (B.20) 

c g B a " 

g4=t 

The correlator is a contraction of X^'^^ and SvY^'^^ summed only over the time shces 
t and t + 1. 

Note that the field vectors X^"^) and SvY^'^^ are independent of t and depend only 
on the distance vector of the bilinear, V = S + J-'. Thus, once they are known, we can 
calculate all bilinear correlators of distance vector V at every time separation with 
no additional inversions. 
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APPENDIX C 



SU{3) PROJECTION 



The application of hypercubic blocking, as discussed in Section [5.6| , requires the 
projection of an arbitrary 3x3 matrix onto the group SU{3). We describe here the 
algorithm used for this projection. For clarity we will use capital letters to denote 
3x3 matrices and lower-case letters to denote 2x2 matrices. 

Given the 3x3 complex matrix M, we wish to find the nearest SU (3) matrix G: 

proj [M] ^ G (C.l) 

SU{3) 

We define the nearest group element to be the one which maximizes: 

tr[GMt] (C.2) 

In order to maximize this trace, we must choose a G such that GM^" is as nearly pro- 
portional to identity as possible. For SU (2) this problem can be solved in closed form. 
While this is not true for SU{3), we can break the problem down into an iterative 
procedure of repeatedly applying the closed-form solution on SU{2) subgroups. If we 
use a set of subgroups which span the full SU{3) group, the process will converge on 
the correct G. 

We begin with a guess for G: 

Gi = l (C.3) 
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This then gives us an initial residual matrix R: 

Ri = GiM^ = 



(C.4) 



The residual matrix is the matrix we wish to make proportional to identity. 

We now enter an iterative process in which each iteration pushes R closer to 
identity and leads us to a more refined value for G. At the start of each iteration we 
choose an SU{2) subgroup to work within. In practice we choose from among three 
SU (2) subgroups, using each in turn. We extract the appropriate 2x2 matrix r from 
R: 



r — r 



(a) 



.a) 



Rn;ll 


Rn;12 


r(2) = 


Rn;ll 


Rn;13 


^(3) ^ 


Rn;22 


Rn;23 


_Rn;21 


Rn;22_ 




Rn;31 


Rn;33 




_Rn;32 


Rn;33_ 



(C.5) 
(C.6) 



where a identifies the subgroup chosen. 

In the context of the SU (2) subgroup, there exists a closed-form expression for 
the group element u nearest the matrix r. We calculate directly the unnormalized 



coefficients o;^ of the matrix u: 



1. 



(54 = Re<j -Tr[r] 



where a, are the Pauli matrices: 



Re<j --Tr[rai] 



(C.7) 



1 

1 



cr2 



-i 

1 



1 
-1 



(C.8) 



After normalizing the coefficients: 



a. 



a. 



\a\ 



u is constructed: 



u 



(C.9) 



(C.IO) 
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Clearly, if u is the closest SU{2) group element to r, applying the inverse of to r 
will result in a matrix as near to proportional to identity as is possible. 

We now return to the full 3x3 matrices in order to complete the iteration: 

a = 1 



U 



( 


Mil 


U12 u 






U22 




_ 


1 






Ui2 


< 





1 






1*22 




"1 


" 







Uu Ui2 


< 





U21 U22_ 



(C.ll) 



a — 3 

We move R closer to identity, and refine our value for G, by applying to G the inverse 
of U. As U is unitary, its inverse is simply its adjoint: 



Gn+l — U^Gn 



(C.12) 
(C.13) 



At this point we have completed one iteration. The process is now repeated using 
our refined guess for G. 

When the trace of R stops improving, we have reached our final value for G. We 
end the iterative process when: 



Tr Rn+i — Tr it!„ < e 



(C.14) 



where e is our error tolerance. In practice we only perform this test after using all 
three SU (2) subgroups in turn, as it is possible that the trace may be flat with respect 
to one subgroup, while improvement is still possible along the other directions. 
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APPENDIX D 
FIGURES 




Figure D.l: Markov chain of ensemble A with thermahzation point Nt and block length 
Nb shown. The correlator C^^^-t is calculated using my = 0.01. Of the ensemble's two 
Markov chains, the chain which begins above the equilibrium value has a disordered initial 
condition. The second has an ordered initial condition. 
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Figure D.2: Markov chain of ensemble B with thermalization point Nt and block length 
Nb shown. The correlator C^^^-t is calculated using niv = 0.01. The Markov chain has an 
ordered initial condition. 
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Figure D.3: Markov chain of ensemble C with thermalization point Nx and block length 
Nb shown. The correlator C^^^-t is calculated using my = 0.01. The Markov chain has an 
ordered initial condition. 
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Figure D.4: Markov chain of ensemble W with thermalization point Nt and block length 
Nb shown. The correlator C^^^-t is calculated using niv = 0.01. The Markov chain has an 
ordered initial condition. 
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Figure D.5: Markov chain of ensemble X with thermahzation point Nt and block length 
Nb shown. The correlator C^^^-t is calculated using my = 0.01. The Markov chain has an 
ordered initial condition. 
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Figure D.6: Markov chain of ensemble Y with thermalization point Nt and block length 
Nb shown. The correlator C^^^-t is calculated using my = 0.01. The Markov chain has an 
ordered initial condition. 
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Figure D.7: Markov chain of ensemble Z with thermalization point Nt and block length 
Nb shown. The correlator C^^^-t is calculated using my = 0.01. The Markov chain has an 
ordered initial condition. 
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Figure D.8: Effective potential for ensembles A hyp and A. The minimum time separation 
chosen is tmin = 4. 
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Figure D.9: Effective potential for ensembles B hyp and B. The minimum time separation 
chosen is tmin = 4. 
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Figure D.IO: Effective potential for ensembles C hyp and C. The minimum time separation 
chosen is tmin = 4. 
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Figure D.ll: Effective potential for ensembles W hyp and W. The minimum time separa- 
tion chosen is tmin = 2. 
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Figure D.12: Effective potential for ensembles X hyp and X. The minimum time separation 
chosen is tmin = 2. 
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Figure D.13: Effective potential for ensembles Y hyp and Y. The minimum time separation 
chosen is tmin = 2. 
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Figure D.14: Effective potential for ensembles Z hyp and Z. The minimum time separation 



chosen is imin = 2. 
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Figure D.15: Effective potential for ensemble Q hyp. The minimum time separation chosen 

is tmin — 4. 
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Figure D.16: Static quark potential for ensemble A hyp. The x's correspond to the un- 
corrected static quark potential, while the diamonds correspond to the corrected potential. 
The result is ro/a = 3.570(27). 
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Figure D.17: Static quark potential for ensemble B hyp. The x's correspond to the un- 
corrected static quark potential, while the diamonds correspond to the corrected potential. 
The result is ro/a = 3.684(44). 
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Figure D.18: Static quark potential for ensemble C hyp. The x's correspond to the un- 
corrected static quark potential, while the diamonds correspond to the corrected potential. 
The result is ro/a = 4.248(78). 
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Figure D.19: Static quark potential for ensemble A. The x 's correspond to the uncorrected 
static quark potential, while the diamonds correspond to the corrected potential. The result 
is ro/a = 3.9(10). 
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Figure D.20: Static quark potential for ensemble B. The x's correspond to the uncorrected 
static quark potential, while the diamonds correspond to the corrected potential. The result 
is ro/a = 4.1(17). 



229 



c 



1.6 


' 1 ' 


1 


1 


1 ' 


1.2 








_ 


o 




















0.8 










0.4 


/ / , 1 


ro/a 
1 


1 


1 



J I L 

2 4 

S 



Figure D.21: Static quark potential for ensemble C. The x's correspond to the uncorrected 
static quark potential, while the diamonds correspond to the corrected potential. The result 
is ro/a = 4.2(12). 
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Figure D.22: Static quark potential for ensemble W hyp. The x's correspond to the un- 
corrected static quark potential, while the diamonds correspond to the corrected potential. 
The result is ro/a = 1.856(19). 



231 




Figure D.23: Static quark potential for ensemble X hyp. The x's correspond to the un- 
corrected static quark potential, while the diamonds correspond to the corrected potential. 
The result is ro/a = 1.847(17). 
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Figure D.24: Static quark potential for ensemble Y hyp. The x's correspond to the un- 
corrected static quark potential, while the diamonds correspond to the corrected potential. 
The result is ro/a = 1.885(19). 
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Figure D.25: Static quark potential for ensemble Z hyp. The x's correspond to the un- 
corrected static quark potential, while the diamonds correspond to the corrected potential. 
The result is ro/a = 1.905(19). 
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Figure D.26: Static quark potential for ensemble W. The x 's correspond to the uncorrected 
static quark potential, while the diamonds correspond to the corrected potential. The result 
is ro/a = 2.10(30). 
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Figure D.27: Static quark potential for ensemble X. The x's correspond to the uncorrected 
static quark potential, while the diamonds correspond to the corrected potential. The result 
is ro/a = 1.78(21). 
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Figure D.28: Static quark potential for ensemble Y. The x's correspond to the uncorrected 
static quark potential, while the diamonds correspond to the corrected potential. The result 
is ro/a = 1.38(17). 
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Figure D.29: Static quark potential for ensemble Z. The x's correspond to the uncorrected 
static quark potential, while the diamonds correspond to the corrected potential. The result 
is ro/a = 1.51(19). 
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Figure D.30: Static quark potential for ensemble Q hyp. The x's correspond to the un- 
corrected static quark potential, while the diamonds correspond to the corrected potential. 
The result is ro/a = 3.650(19). 
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Figure D.31: Dependence of the Sommer scale on tmin for ensembles A hyp, A, B hyp, B, 
C hyp, and C. The filled diamond corresponds to tmm used in the final fit. 
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Figure D.32: Dependence of the Sommer scale on tmm for ensembles W hyp, W, X hyp, 
and X. The filled diamond corresponds to imin used in the final fit. 
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Figure D.33: Dependence of the Sommer scale on tmin for ensembles Y hyp, Y, Z hyp, and 
Z. The filled diamond corresponds to tmin used in the final fit. 
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Figure D.34: Dependence of the Sommer scale on fmin for ensemble Q hyp. The filled 
diamond corresponds to imin used in the final fit. 
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Figure D.35: Dependence of the Sommer scale on tmax for ensembles A hyp, A, B hyp, B, 
C hyp, and C. The filled diamond corresponds to tmax used in the final fit. 
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Figure D.36: Dependence of the Sommer scale on tmax for ensembles W hyp, W, X hyp, 
and X. The filled diamond corresponds to tmax used in the final fit. 
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Figure D.37: Dependence of the Sommer scale on tmax for ensembles Y hyp, Y, Z hyp, and 
Z. The filled diamond corresponds to tmax used in the final fit. 
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Figure D.38: Dependence of the Sommer scale on tmax for ensemble Q hyp. The filled 
diamond corresponds to imax used in the final fit. 
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Fi gure D.39: Dependence of the Sommer scale on Smin for ensembles A hyp, A, B hyp, B, 
C hyp, and C. The filled diamond corresponds to Smax used in the final fit. 
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Figure D.40: Dependence of the Sommer scale on Smin for ensembles W hyp, W, X hyp, 
and X. The filled diamond corresponds to Smax used in the final fit. 
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Figure D.41: Dependence of the Sommer scale on -Smin ensenibles V hyp, Y, Z hyp, and 
Z. The filled diamond corresponds to Smax used in the final fit. 
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Figure D.42: Dependence of the Sommer scale on Smin for ensemble Q hyp. The filled 
diamond corresponds to Smax used in the final fit. 
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Figure D.43: Effective pion mass for ensembles A hyp, A, B hyp, B, C hyp, and C at 
my = 0.01. The minimum time separation chosen is tmin = 7. 
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Figure D.44: Effective pion mass for ensembles W hyp, W, Z hyp, and Z at mv = 0.01. 
The minimum time separation chosen is tmin = 4. 
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Figure D.45: Effective pion mass for ensembles Y hyp, Y, Z hyp, and Z at my = 0.01. The 
minimum time separation chosen is imin = 4. 
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Figure D.46: Effective pion mass for ensemble Q hyp at my = 0.01. The minimum time 
separation chosen is tmin = 7. 
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Figure D.47: Quadratic pion-mass fit and kaon quark-mass threshold for ensembles A hyp, 
A, B hyp, B, C hyp, and C. The dotted vertical lines correspond to the determined values 
of niQj,. 
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Figure D.48: Quadratic pion-mass fit and kaon quark-mass tlireshold for ensembles W hyp, 
W, X hyp, and X. Tlie dotted vertical lines correspond to the determined values of mq^^ . 
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Figure D.49: Quadratic pion-mass fit and kaon quark-mass thresliold for ensembles Y hyp, 
Y, Z hyp, and Z. The dotted vertical lines correspond to the determined values of rnQ^^. 
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Figure D.50: Quadratic pion-mass fit and kaon quark-mass threshold for ensemble Q hyp. 
The dotted vertical line corresponds to the determined value of mq^ . 
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Figure D.51: P^r degree of freedom over a range of valence-quark-mass cutoffs for 
ensemble A hyp. The filled diamond corresponds to the cutoff used in the final fit. 
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Figure D.52: Pion mass squared versus valence quark mass from the pqChPT fit of ensem- 
ble A hyp. Filled diamonds correspond to valence-quark-mass values within the fit range, 
while open diamonds correspond to those values beyond it. 
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Figure D.53: Pion decay constant versus valence quark mass from the pqChPT fit of 
ensemble A hyp. Filled diamonds correspond to valence-quark-mass values within the fit 
range, while open diamonds correspond to those values beyond it. 
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Figure D.54: Rm versus valence quark mass from the pqCliPT fit of the correlator of 
ensemble A hyp. Filled diamonds correspond to valence-quark-mass values within the fit 
range, while open diamonds correspond to those values beyond it. 
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Figure D.55: Results from the pqChPT fit of ensemble A hyp. Filled diamonds correspond 
to valence-quark-mass values within the fit range, while open diamonds correspond to those 
values beyond it. 
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Figure D.56: Results from the pqChPT fit of ensemble B hyp. Filled diamonds correspond 
to valence-quark-mass values within the fit range, while open diamonds correspond to those 
values beyond it. 
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Figure D.57: Results from the pqChPT fit of ensemble A. Filled diamonds correspond to 
valence-quark-mass values within the fit range, while open diamonds correspond to those 
values beyond it. 
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Figure D.58: Results from the pqChPT fit of ensemble B. Filled diamonds correspond to 
valence-quark-mass values within the fit range, while open diamonds correspond to those 
values beyond it. 
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Figure D.59: Results from the pqChPT fit of ensemble W hyp. Filled diamonds correspond 
to valence-quark-mass values within the fit range, while open diamonds correspond to those 
values beyond it. 
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Figure D.60: Results from the pqChPT fit of ensemble X hyp. Filled diamonds correspond 
to valence-quark-mass values within the fit range, while open diamonds correspond to those 
values beyond it. 
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Figure D.61: Results from the pqChPT fit of ensemble Y hyp. Filled diamonds correspond 
to valence-quark-mass values within the fit range, while open diamonds correspond to those 
values beyond it. 
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Figure D.62: Results from the pqChPT fit of ensemble Z hyp. Filled diamonds correspond 
to valence-quark-mass values within the fit range, while open diamonds correspond to those 
values beyond it. 
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Figure D.63: Results from the pqChPT fit of ensemble W. Filled diamonds correspond to 
valence-quark-mass values within the fit range, while open diamonds correspond to those 
values beyond it. 
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Figure D.64: Results from the pqChPT fit of ensemble X. Filled diamonds correspond to 
valence-quark-mass values within the fit range, while open diamonds correspond to those 
values beyond it. 
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Figure D.65: Results from the pqChPT fit of ensemble Y. Filled diamonds correspond to 
valence-quark-mass values within the fit range, while open diamonds correspond to those 
values beyond it. 
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Figure D.66: Results from the pqChPT fit of ensemble Z. Filled diamonds correspond to 
valence-quark-mass values within the fit range, while open diamonds correspond to those 
values beyond it. 
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Figure D.67: Results from the pqChPT fit of ensemble Q hyp. Filled diamonds correspond 
to valence-quark-mass values within the fit range, while open diamonds correspond to those 
values beyond it. 
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Figure D.69: Rm versus valence quark mass from the cubic fit of the correlator of ensemble 
A hyp. 
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Figure D.70: Results from the pqChPT fit of ensemble C hyp. Filled diamonds correspond 
to valence-quark-mass values within the fit range, while open diamonds correspond to those 
values beyond it. 
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Figure D.71: Results from the pqChPT fit of ensemble C. Filled diamonds correspond to 
valence-quark-mass values within the fit range, while open diamonds correspond to those 
values beyond it. 
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Figure D.72: Dependence of the results of the pqChPT fit of ensemble A hyp on the 
valence-quark-mass cutoff. The filled diamond corresponds to the cutoff used in the final 
fit. 
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Figure D.73: Dependence of the results of the pqChPT fit of ensembles A hyp, A, B hyp, 
B, C hyp, and C on the valence-quark-mass cutoff. The filled diamond corresponds to the 
cutoff used in the final fit. 
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Figure D.74: Dependence of the results of the pqChPT fit of ensembles W hyp, Z, W hyp, 
and W on the valence-quark-mass cutoff. The filled diamond corresponds to the cutofi^ used 
in the final fit. 
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Figure D.75: Dependence of the results of the pqChPT fit of ensembles Y hyp, Y, Z hyp, 
and Z on the valence-quark-mass cutoff. The filled diamond corresponds to the cutoff used 
in the final fit. 
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Figure D.76: Dependence of the results of the pqChPT fit of ensemble Q hyp on the 
valence-quark-mass cutoff. The filled diamond corresponds to the cutoff used in the final 
fit. 
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Figure D.77: Rm versus valence quark mass from the simultaneous pqChPT fit of the 
correlators of ensembles W hyp through Z hyp. Filled diamonds correspond to valence- 
quark-mass values within the fit range, while open diamonds correspond to those values 
beyond it. An average of the ensembles' values for mg^ is used. 
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Figure D.78: Pion decay constant versus valence quark mass from the simultaneous 
pqChPT fit of ensembles W hyp through Z hyp. Filled diamonds correspond to valence- 
quark-mass values within the fit range, while open diamonds correspond to those values 
beyond it. An average of the ensembles' values for mg^ and lattice spacing are used. 
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Figure D.79: Pion decay constant versus valence quark mass from the simultaneous 
pqChPT fit of ensembles W through Z. Filled diamonds correspond to valence-quark-mass 

values within the fit range, while open diamonds correspond to those values beyond it. An 
average of the ensembles' values for mg^ and lattice spacing are used. 
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Figure D.80: Pion decay constant versus valence quark mass from the simultaneous 
pqChPT fit of ensembles W through Z. Filled diamonds correspond to valence-quark-mass 

values within the fit range, while open diamonds correspond to those values beyond it. An 
average of the ensembles' values for mg^ and lattice spacing are used. 
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Figure D.81: Rm versus valence quark mass from the four independent pqCliPT fits of 
ensembles W hyp through Z hyp. Filled diamonds correspond to valence-quark-mass values 
within the fit range, while open diamonds correspond to those values beyond it. An average 
of the ensembles' values for mg^ is used. 
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Figure D.82: Pion decay constant versus valence quark mass from the four independent 
pqChPT fits of ensembles W hyp through Z hyp. Filled diamonds correspond to valence- 
quark-mass values within the fit range, while open diamonds correspond to those values 
beyond it. An average of the ensembles' values for mg^ and lattice spacing are used. 
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Figure D.83: Rm versus valence quark mass from the four independent pqChPT fits of 
ensembles W through Z. Filled diamonds correspond to valence-quark-mass values within 
the fit range, while open diamonds correspond to those values beyond it. An average of the 
ensembles' values for mg^ is used. 
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Figure D.84: Pion decay constant versus valence quark mass from the four independent 
pqChPT fits of ensembles W through Z. Filled diamonds correspond to valence-quark-mass 

values within the fit range, while open diamonds correspond to those values beyond it. An 
average of the ensembles' values for mg^ and lattice spacing are used. 
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Figure D.85: Results from the simultaneous pqChPT fit of ensembles W hyp through Z hyp, 
displayed along the unquenched line of the quark-mass plane. An average of the ensembles' 
values for mg^ and lattice spacing are used. 
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Figure D.86: Results from the simultaneous pqCliPT fit of ensembles W through Z, dis- 
played along the unquenched line of the quark-mass plane. An average of the ensembles' 
values for mg^ and lattice spacing are used. 
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Figure D.87: Dependence of the results of the simultaneous pqChPT fit of ensembles 
W hyp through Z hyp on the valence-quark-mass cutoff. The filled diamond corresponds to 
the cutoff used in the final fit. An average of the ensembles' values for is used. 
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Figure D.88: Dependence of the results of the simultaneous pqChPT fit of ensembles W 
through Z on the valence-quark-mass cutoff. The filled diamond corresponds to the cutoff 
used in the final fit. An average of the ensembles' values for ttiq^ is used. 
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